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ABSTRACT 

We present results from thirteen cosmological simulations that explore the parameter space 
of the “Evolution and Assembly of GaLaxies and their Environments” (EAGLE) simulation 
project. Eour of the simulations follow the evolution of a periodic cube L = 50 cMpc on a 
side, and each employs a different subgrid model of the energetic feedback associated with 
star formation. The relevant parameters were adjusted so that the simulations each reproduce 
the observed galaxy stellar mass function at z = 0.1. Three of the simulations fail to form 
disc galaxies as extended as observed, and we show analytically that this is a consequence of 
numerical radiative losses that reduce the efficiency of stellar feedback in high-density gas. 
Such losses are greatly reduced in the fourth simulation - the EAGLE reference model - by 
injecting more energy in higher density gas. This model produces galaxies with the observed 
size distribution, and also reproduces many galaxy scaling relations. In the remaining nine 
simulations, a single parameter or process of the reference model was varied at a time. We 
hnd that the properties of galaxies with stellar mass < M* (the “knee” of the galaxy stellar 
mass function) are largely governed by feedback associated with star formation, while those of 
more massive galaxies are also controlled by feedback from accretion onto their central black 
holes. Both processes must be efficient in order to reproduce the observed galaxy population. 
In general, simulations that have been calibrated to reproduce the low-redshift galaxy stellar 
mass function will still not form realistic galaxies, but the additional requirement that galaxy 
sizes be acceptable leads to agreement with a large range of observables. 
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1 INTRODUCTION 

The formation, assembly and evolution of cosmic structures is or¬ 
chestrated by gravitational collapse. The non-linearity of this pro¬ 
cess precludes a fully-predictive analytic theory of structure for¬ 
mation, requiring that the confrontation of theoretical expectations 
with observational measurements must generally proceed via nu¬ 
merical simulations. The predictions of cosmological simulations 
based on the prevailing A-cold dark matter (ACDM) paradigm, in 


the regime where those outcomes are determined primarily by grav¬ 
itational forces, have been corroborated by a diverse range of ob¬ 
servational tests. These include, but are not limited to, cosmic shear 
induced by large-scale structure (e.g. |Fu e t al.|2008 ), the ab undance 
of brightest cluster galaxies (BCGs, e.g. Rozo et al.|2010^ , tests of 
the cosmic expansion rate (e.g. Blake etal.|201 la l and the distance- 
redshift relation (e.g. |Blake et ^ 201 Ib^ , redshift-space distortions 
of the 2-point correlation function (e.g. |Beutler et al.||201^ and 
the luminosity-distance relation of type la supernovae (e.g. Suzukij 
let al.|2012| >. 


* E-mail: r.a.crain@ljmu.ac.uk 


The formation and evolution of galaxies is governed ulti- 
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mately, however, by the interaction of the diverse physical pro¬ 
cesses that, in addition to gravity, influence baryonic matter. The 
inclusion of these processes in simulations is recognised as a major 
challenge, owing both to the complexity of the physical processes, 
and the difficulty of developing numerical algorithms able to ac¬ 
curately model their effects in a computationally efficient manner. 
This challenge has, by and large, impeded cosmological hydrody- 
namical simulations from yielding galaxy populations whose prop¬ 
erties are consistent with observational measurements. Although 
imperfect models can prove instructive, greater confidence is gen¬ 
erally ascribed to those that more accurately resemble the observed 
Universe. Moreover, the reproduction of key observables is often a 
prerequisite for testing particular aspects of galaxy formation the¬ 
ory. For example, when wishing to study the evolution of angular 
momentum in disc galaxies, a model that reproduces their observed 
size and rotation velocity is clearly desirable. 


The reproduction of these particular diagnostics has in fact 
become a cause celebre for the simulation community, owing to 
the long-standing need to address the closely related “overcooling” 


(Cole|1991||White & Frenk|1991[|Blanchard et al.|1992 

IBalogh 

et al.|2001b and “angular momentum” (|Katz & Gunn|1991| 

Navarro 

& White 1994^ problems. In the absence of feedback. 

gas effi- 


ciently radiates the heat it acquires from thermalising its gravita¬ 
tional potential. This excess dissipation has two principal conse¬ 
quences: i) the fraction of gas that is converted into stars by the 
present epoch is much higher than observed, and ii) the forma¬ 
tion and coalescence of dense clumps spuriously drains angular 
momentum from the baryons. Simulated galaxies therefore form 
too many stars (and do so too early), they are more compact than 
observed, and they exhibit insufficient rotational support. The in¬ 
clusion of prescriptions for energetic feedback processes in mod¬ 
els has been shown to allevi ate these problems jA badi et a l.|2003[ 
|Sommer-Larsen et al.|2003l|Springel & Hernquist|2003| , and has 

enabled several groups to conduct simulations of the ACDM cos¬ 
mogony that form galaxies with sizes and rotation curves that are, 
for particular galaxy masses, consistent with observational mea- 


surements (e.g. Govemato et al.|2004 |Okamoto et al. 

20051 Sales 

et al. 

2010| Guedes et al. 

201 1| McCarthy et al.|2012| 

Brook et al. 

2012 

|Munshietal.|20131 

Aumer et al.|2013[|Marinacci et al.|2014b. 


In spite of this success, the detailed behaviour of the multi¬ 
phase interstellar medium (ISM) when subject to energetic feed¬ 
back remains ill-understood, and the community has yet to con¬ 
verge on unique solutions to the overcooling and angular momen¬ 
tum problems ( [Scannapieco et al.|2012^ . The principal uncertainty 
is arguably one of accounting. Firstly, it is not known what are the 
energy, momentum and mass fluxes incident upon the ISM and star¬ 
forming complexes therein (but see|Lopez et al.|201 l||Rosen et al.| 
|2014[ l, due to mechanisms such as radiation pressure and winds 
from O-class stars, asymptotic giant branch (AGB) stars and active 
galactic nuclei (AGN); the photoionisation and photoelectric heat¬ 
ing of HII regions by radiation associated with stars (including X- 
ray binaries) and the accretion discs of black holes (BHs); and ther¬ 
monuclear and core collapse supemovae (SNe). A second, often 
overlooked issue is that it is unknown what fraction of the incident 
energy is dissipated by radiative processes and thermal conduction 
(e.g. |Orlando et al.|2005^ , and what fraction of the incident momen¬ 
tum is lost due to cancellation. Estimating these initial “losses” is a 
long-standing problem in the study of the ISM, not least because of 
the extreme resolution and dynamic range demands of the problem: 
the losses are typically established on scales significantly smaller 
than a parsec (e.g. |Me l lema et al.|2002[ [Fragile et al.|200 4[ [Yirak| 
|et aLpOfO]). This is at least three orders of magnitude smaller than 


the typical size of an ISM resolution element in simulations of large 
cosmological volumes. 

Since these losses cannot be modelled directly by cosmolog¬ 
ical simulations, their impact on resolved scales must be incorpo¬ 
rated into phenomenological “subgrid” treatments that approximate 
the action of unresolved processes, and couple them to resolved 
scaleQ The implementation and parametrisation of subgrid rou¬ 
tines is therefore the greatest source of uncertainty in cosmological 
simulations, and adjustment of these characteristics can result in the 


Schaye et al.|2010 

|Haas et al.|2013a|b 

[Scannapieco et al. 

2012 

Vogelsberger et al. 

2013||Le Brun et al. 

2014[|Torrey et al.|2014|( 


Until small-scale losses can be accurately computed and appropri¬ 
ately incorporated into subgrid routines, it will remain impossible 
to formulate a truly predictive cosmological simulation that can, for 
example, yield ab initio estimates of the stellar mass of galaxies or 
the mass of their central BH. 

In a companion paper, [Schaye et aL] ( |20I5| hereafter SI5) ar¬ 
gue that the appropriate methodology is therefore to calibrate the 
parameters of subgrid routines, in order that simulations reproduce 
well characterised observables. The calibrated observables them¬ 
selves cannot then be advanced as predictions of the model, but 
those not considered during the calibration can reasonably be con¬ 
sidered as consequences of the implemented astrophysics. An ob¬ 
vious advantage of this approach is that, by ensuring that key prop¬ 
erties of the galaxy population are reproduced, simulations can be 
used to address the widest range of problems. A related advantage 
is that, since any alteration to the resolution of a calculation will in 
general necessitate a recalibration of the model, the adopted sub¬ 
grid routines need not sacrifice physical detail in order to realise 
numerical convergence. 

Adopting this philosophy, S15 introduced the “Evolution and 
Assembly of GaLaxies and their Environments” (EAGLE) project, 
a suite of cosmological hydrodynamical simulations of the ACDM 
cosmogony conducted by the Virgo Consortiunj^ Feedback from 
star formation and AGN is implemented thermally, such that out¬ 
flows develop as a result of pressure gradients and without the need 
to impose winds ‘by hand’, for example by specifying their veloc¬ 
ity and mass loading with respect to the star formation rate. The 
parameters of the subgrid routines governing feedback associated 
with star formation and the growth of BHs are calibrated to repro¬ 
duce the observed z = 0.1 galaxy stellar mass function (GSMF) 
and the relation between the mass of galaxies and their central 
BH, respectively, whilst also seeking to yield galaxies with sizes 
(i.e. effective radius) similar to those observed. S15 focussed on 
the EAGLE reference model (“Ref”), and a complementary model 
designed to meet the calibration criteria at higher resolution (“Re- 
caF’jH Besides demonstrating that cosmological simulations can 
be calibrated to reproduce these diagnostics successfully with un¬ 
precedented accuracy, the study showed that the simulations repro¬ 
duce a diverse and representative set of low-redshift observables 
that were not considered during the calibration process. In a sepa¬ 
rate paper, [Furlong et ar](|20I4^ show that the EAGLE simulations 


^ We refer to losses on these scales as “subgrid losses”. Losses induced 
by processes acting on scales that are resolved by cosmological simulations 
can also be significant, and dependent upon the subgrid implementation; we 
term these “macroscopic losses”. 

^ See also http://eagle.strw.leidenuniv.nl and http://icc.dur.ac.uk/Eagle 
® SIS also introduced a third model that better reproduces the observed 
properties of intragroup gas at intermediate resolution by adopting a higher 
AGN heating temperature (“AGNdT9”). 
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also broadly reproduce the observed GSMF as early as 2 = 7, and 
accurately track its evolution to the present day. 

This paper introduces many more simulations from the EA¬ 
GLE project. The key aims of this study are to illustrate the reasons 
for the parametrisation adopted by the EAGLE reference model, 
and the sensitivity of its outcomes to the variation of the key sub¬ 
grid parameters. The simulations explored here are naturally di¬ 
vided into two categories. The first comprises four simulations cal¬ 
ibrated to yield the 2 = 0.1 GSMF and central BIT masses as a 
function of galaxy stellar mass. The models, one of which is the 
EAGLE reference model, differ in terms of the adopted subgrid 
efficiency of feedback associated with star formation, and the fash¬ 
ion by which this efficiency depends (if at all) upon the properties 
of the local environment. The successful reproduction of the cali¬ 
bration diagnostics by each of the models highlights that these ob¬ 
servables alone do not identify a unique “solution”, and indicates 
that complementary constraints are necessary to break modelling 
degeneracies and, potentially, motivate the inclusion of additional 
complexity. 

The second category comprises simulations each featuring a 
variation of a single subgrid parameter value with respect to the 
reference model. These calculations enable an examination of the 
role of these parameters in a fashion similar to the OWLS project 
( [Schaye et aL||201^ , and highlight the sensitivity of outcomes to 
the variation of these parameters. In common with complementary 
studies, these simulations indicate that the properties of the simu¬ 
lated galaxy population are most sensitive to the subgrid parameters 
governing the efficiency of energy feedback (|Schaye^^e£^ |20M 


Scannapieco et al.||2012| |Haas et al.||2013a|b| Vogelsberger et al.| 

2013| l. 

This paper is structured as follows. The simulation initial con¬ 
ditions, and the algorithms used to evolve them, are described in 
The parametrisation of the four models that are calibrated to re¬ 
produce the 2 = 0.1 GSME is described in Results from these 
simulations, which serve as a motivation for the development of 
the reference model, are presented in ^ where results from simu¬ 
lations featuring single-parameter variations of the reference model 
are also shown. Finally, the results are summarised and discussed 


2 SIMULATIONS AND SUBGRID PHYSICS 

This section comprises an overview of the simulation setup and 
subgrid physics implementation. It includes similar information to 
Sections 3 & 4 of SI5, so readers familiar with the simulations may 
wish to skip this section. A relatively comprehensive description 
of the subgrid implementations of star formation and feedback is 
retained here, because these details are a necessary foundation for 
later sections. 

The cosmological parameters assumed by the EAGLE sim¬ 
ulations are those recently inferred by the jPlanck Collaboration! 
( |2014a|b| >, the key parameters being fim = 0.307, Ha = 0.693, 
fib = 0.04825, h = 0.6777 and a% = 0.8288. Initial conditions 
adopting these parameters were g enerated using a tr ansfer function 
created with the CAMB software jLewis et al.|2000) , the 2"‘^-order 
Lagrangian perturbation theory method of[Jenkii^^20]T)|, and the 
Gaussian white noise field Panphasia ([Jenkins |2013[ Jenkins &| 
|Booth|2013[. A complete description of the generation of the initial 


conditions is provided in Appendix B of S15, and the tools neces¬ 
sary to generate them independently are available onlin^ 

The simulations were evolved by a modified version of the 
A-body TreePM smoothed particle hydrodynamics (SPH) code 
Gadget3, last described by [Spring^ j2005) . The modifications 
comprise updates to the hydrodynamics algorithm and the time¬ 
stepping criteria, and the addition of subgrid routines governing 
the phenomenological implementation of processes occurring on 
scales below the resolution limit of the simulations. The updates 
to the hydrodynamics algorithm, which we collectively refer to as 
“Anarchy” (Dalla Vecchia in prep.), comprise an implementation 
of the pressure-entropy formulation of SPH derived by [Hopkins 


(2013 


(2010 


Price 


the artificial viscosity switch proposed by [Cullen & Dehnm 
an artificial conduction switch similar to that proposed by 


2008 1 , the C smoothing kernel of 


time-step limiter of[Durier & Dalla Vecc 


Wendland 

< 2012 ) . 


11995 I, and the 


The subgrid routines represent an evolution of those used for 
the GIMIC ( [Crain et al.[|200^, OW LS jSchaye et ai:|[20T()| l and 
cosmo-OWLS <Le Brun et al.[2014) projects, and include element- 
by-element radiative cooling and photoionisation heating for 11 
species, star formation, stellar mass loss, energy feedback from star 
formation, gas accretion onto and mergers of BHs, and AGN feed¬ 
back. The key updates with respect to the routines used by OWLS 
are the inclusion of a metallicity dependence in the star formation 
law, the implementation of energy feedback associated with star 
formation via stochastic thermal heating, and the inclusion of a vis¬ 
cous transport limit on the BH accretion rate. 

S15 introduced the resolution nomenclature of the EAGLE 
project. “Intermediate-resolution” simulations have particle masses 
corresponding to an L = 100 comoving Mpc (hereafter cMpc) 
volume realised with 2 x 1504® particles (an equal number of 
baryonic and dark matter particles), such that the initial gas par¬ 
ticle mass is mg = 1.81 x 10® Mq, and the mass of dark mat¬ 
ter particles is rudm = 9.70 x 10® M©. The Plummer-equivalent 
gravitational softening length is fixed in comoving units to 1/25 
of the mean interparticle separation (2.66 comoving kpc, here¬ 
after ckpc) until 2 = 2.8, and in proper units (0.70 proper kpc, 
hereafter pkpc) at later times. The intermediate-resolution simu¬ 
lations marginally resolve the Jeans scales at the star formation 
threshold (nn 22 10“^ cm“®) in the warm (T ~ lO'^K) ISM. 
“High-resolution” simulations adopt particle masses and softening 
lengths that are smaller by factors of eight and two, respectively. 
The SPH kernel size, specifically its support radius, is limited to 
a minimum of one-tenth of the gravitational softening scale. This 
study focusses on intermediate-resolution simulations using vol¬ 
umes of side L = 25, 50 and 100 cMpc, which therefore comprise 
2 X 376®, 2 X 752® and 2 x 1504® particles, respectively. 

Galaxies and their host haloes are identified by a multi-stage 
process, beginning with the application of the friends-of-friends 
(EoF) algorithm ( [Davis et al.[1985) to the dark matter particle distri¬ 
bution, with a linking length of 6 = 0.2 times the mean interparticle 
separation. Gas, star and BH particles are associated with the EoE 
group, if any, of their nearest neighbour dark matter particle. The 
SUBFIND algorithm ( [Springel et aTjllOOT] [Dolag et al.[[200^ is 
then used to identify self-bound substructures, or subhaloes, within 
the full particle distribution (gas, stars, BHs and dark matter) of 
FoF haloes. The subhalo comprising the particle with the minimum 
gravitational potential, which is almost exclusively the most mas¬ 
sive subhalo, is defined as the central subhalo, the remainder being 


See http://eagle.strw.leidenuniv.nl 
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satellite subhaloes. The coordinate of the particle with the mini¬ 
mum potential also defines the position of the halo, about which 
is computed the spherical overdensity (SO; [Lacey & Co!^|1994^ 
mass, M 20 Q, for the adopted density contrast of 200 times the crit¬ 
ical density, pc- Satellite subhaloes separated from their central 
galaxy by less than the minimum of 3 pkpc and the stellar half¬ 
mass radius of the central galaxy are merged into the latter; this 
step eliminates a small number of low-mass subhaloes dominated 
by single, high-density gas particles or BHs. Finally, when quoting 
the properties of galaxies (e.g. stellar mass, star formation rate), 
only those subhalo particles within a spherical aperture of radius 
30 pkpc are considered. S15 (their Figure 6 ) demonstrated that this 
practice yields stellar masses comparable to those recovered within 
a projected circular aperture with the Petrosian radius at z = 0.1. 


2.1 Radiative processes 


Radiative cooling and heating rates are computed on an element- 
by-element basis by interpolating tables, generated with CLOUDY 
(version 07.02, [Perland et aL]|1998^ , that specify cooling rates as 
a function of density, temperature and redshift, under the assump¬ 
tion that the gas is optically thin, is in ionisation equilibrium, and 
is exposed to the cosmic microwave background and a spatially- 
uniform, temporally-evolving [Flaardt & Madau ( 200 1|> UV/ X-ray 
background (for further details, see Wiersma et al. 2009a^ . The 
UV/X-ray background is imposed instantaneously at z = 11.5. To 
account for enhanced photoheating rates (relative to the optically 
thin rates assumed here) during the epochs of reionisation, 2eV 
per proton mass is injected, rapidly heating gas to ~ 10"^ K. This 
is done instantaneously at z = 11.5 (consistent with Planck con¬ 
straints) for HI reionisation, but for Hell the energy injection is dis¬ 
tributed in redshift with a Gaussian function centred about 2 = 3.5 
with a width of (j{z) = 0.5. This ensures that the thermal evolution 
of the intergalactic medium mimics that inferred by|Schaye et al.| 

( llOOOt . 


2.2 The ISM and star formation 

Simulations of large cosmological volumes lack, in general, the 
resolution and physics to model the cold (T <C 10"^ K) inter¬ 
stellar gas phase from which molecular clouds and stars form. A 
global temperature floor, Teos (p) is therefore imposed, correspond¬ 
ing to a polytropic equation of state, Peos oc p^'“=, normalised 
to Teas = 8000 K at nn = 0.1 cm“®. A fiducial polytrope of 
7 eos = 4/3 is adopted, since this ensures that the Jeans mass, and 
the ratio of the Jeans length to the SPH support radius, are inde¬ 
pendent of density, thus inhibiting spurious fragmentation ( jSchayej 
|& Dalla Vecchia|2008| >. The effect of varying 7 eos is explored in 
§ |4.2| where simulations conducted using isothermal ( 7 eos = 1) 
and adiabatic ( 7 eos = 5/3) equations of state are examined. 

A second temperature floor of 8000 K is imposed for gas with 
riH > 10 “®cm“^, which prevents metal-rich gas from cooling 
to very low temperatures, since the physical processes required to 
model dense, low-temperature gas are not included here. This floor 
does not apply to low-density (i.e. intergalactic) gas, since such gas 
cools adiabatically and is modelled accurately by the hydrodynam¬ 
ics scheme. 

Star formation is implemented stochastically, based on the 
pressure law scheme of jSchaye & Dalla Vecchia| ( |2008^ . Under the 
(reasonable) assumption that star-forming gas is self-gravitating, 
the observed Kennicutt-Schmidt star formation law (|Kennicutt| 


[T^ , 


E* = A 


1 Mq pc“ 


( 1 ) 


where E* and Eg are the surface density density of stars and gas, 
respectively, can be expressed as: 


, _9, —71 / 'y \ t)/^ 

m. = rrigA (1 Mq pc ) > (2) 

where rUg is the gas particle mass, 7 = 5/3 is the ratio of specific 
heats (and should not be confused with 7 eos), G is the gravitational 
constant, /g is the mass fraction in gas (assumed to be unity), and 
P is the total pressure. This pressure law implementation is ad¬ 
vantageous for two reasons. Firstly, the free parameters of the star 
formation law (A, n) are specified explicitly by observations: the 
values A — 1.515 x 10“^ Mq yr“^ kpc“^ and n = 1.4 (n = 2 
for riH > 10® cm“®) are adopted, where the value of A has been 
adjusted from that reported by |Kennicutt| ( | 199^ to convert from the 
Salpeter initial stellar mass function (IMF) to the jChabrierj p003) > 
form adopted by the simulations. Secondly, this implementation 
guarantees that the observed Kennicutt-Schmidt relation is repro¬ 
duced for any equation of state (i.e. any combination of Teas and 
7eos) applied to star-forming gas. This is in contrast to volumet¬ 
ric star formation laws, which must be recalibrated whenever the 
equation of state is altered. 

Star formation occurs only in cold, dense gas, requiring that 
a density threshold for star formation, n^, be imposed. Since the 
transition from a warm, neutral phase to a cold, molecular one 
occurs at lower densities and pressures in metal-rich (and hence 
dust-rich) gas, we adopt the metallicity-dependent star formation 
threshold proposed by jSchayej fl2004|l, which was implemented in 
the OWLS simulation “SFTHRESHZ”: 


n^(Z) = min 


0.1 cm 


0.002 


, 10 cm 


(3) 


where Z is the gas metallicity. Hydrogen number density, nn, is 
related to the overall gas density, p, via nn = Xp/mn, where X 
is the hydrogen mass fraction and mn is the mass of a hydrogen 
atom. To examine the effects, if any, of adopting this metallicity- 
dependent threshold, the EAGLE suite includes a simulation that 
adopts a constant threshold of — 0.1 cm“®, which was the 
fiducial approach of OWLS. In both cases, to prevent star formation 
in low-overdensity gas at high redshift, star-forming gas is required 
to have an overdensity S > 57.7. 


2.3 Stellar evolution and mass loss 

The implementation of stellar evolution and mass loss is based 
upon that described by jWiersma et al.j j2009bj ). Star particles are 
treated as simple stellar populations (SSPs) with an IMF of the form 
proposed by jChabrierj ( j2003T l, spanning the range 0.1 — 100 Mq. 
At each time step and for each stellar particle, those stellar masses 
reaching the end of the main sequence phase are identified using 
metallicity-dependent lifetimes, and the fraction of the initial parti¬ 
cle mass reaching this evolutionary stage is used, together with the 
particle initial elemental abundances and nucleosynthetic yield ta¬ 
bles, to compute the mass of each element that is lost through winds 
from AGB stars, winds from massive stars, and type II SNe. Eleven 
elements are tracked, and the mass and energy lost through type la 
SNe is also computed, assuming the rate of type la SNe per unit 
stellar mass is specified by an empirically-motivated exponential 
delay function. 
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2.4 Energy feedback from star formation 


Stars inject energy and momentum into the ISM via stellar winds, 
radiation, and SNe. These processes are particularly important for 
massive, short-lived stars and, if star formation is sufficiently vig¬ 
orous, the associated feedback can drive large-scale galactic out¬ 
flows (e.g. [Veilleux et al.|[2005l ). At present, simulations of large 
cosmological volumes lack the resolution necessary to model the 
self-consistent development of outflows from feedback injected on 
the scales of individual star clusters, and must appeal to a subgrid 
treatment. 

In the simplest implementation of energy feedback by thermal 
heating, the energy produced at each timestep by a star particle is 
distributed to a number of its neighbouring hydrodynamic resolu- 


tion elements, supplementing their internal energy. Dalla Vecchia 

& Schaye|i 

2012| see also|Dalla Vecchia & Schaye|2008 Creasey 

etal.|20111 

Keller et al.|2014 and Creasey et al.|2015b argue that the 


feedback energy (canonically ~ 10° erg per 100 Mq for a stan¬ 
dard IMF if considering SNe as the sole energy source) is initially 
distributed over too much mass: the mass of at least one resolution 
element, 0(10"^ — 10^) Mq, rather than that of the actual ejecta, 
0{lVp — 10^) Mq. The resulting temperature increment is then far 
smaller than in reality, and by extension the radiative cooling time 
of the heated gas is much too short. Pressure gradients established 
by the heating are too shallow and, perhaps more importantly, are 
typically erased on a (radiative) timescale shorter than the sound 
crossing time of a resolution element. 

The EAGLE simulations adopt the stochastic thermal feed¬ 
back scheme of |Dalla Vecchia & Schaye| ( |2012| l, in which the tem¬ 
perature increment, ATsf, of heated resolution elements is speci¬ 
fied. Besides enabling one to mitigate the problem described above, 
a stochastic implementation of feedback is advantageous because it 
enables the quantity of energy injected per feedback event to be 
specified, even if the mean quantity of energy injected per unit 
stellar mass formed is fixed. Having specified ATsf, the proba¬ 
bility that a resolution element neighbouring a young star particle 
is heated, is determined W the fraction of the energy budget that 
is available for feedbacl^ For consistency with the nomenclature 
introduced by |Dalla Vecchia & Schaye|p012| l, we refer to this frac¬ 
tion as /th. 

We adopt the convention that /th = 1 equates to an ex¬ 
pectation value of the injected energy of 1.736 x 10'*® erg Mg* 
(8.73 X 10*® erg g“*) of stellar mass formed. This corresponds to 
the energy available from type II SNe resulting from a Chabrier 
IMF, subject to two assumptions. Firstly, that 6 — 100 Mq stars are 
the progenitors of type II SNe (6 — SMq stars explode as elec- 
tron capture supemovae in models with convective overshoot; e.g. 


Chiosi et al.|1992 i, and secondly that each SN liberates 10°* erg. 


We inject energy once for each star particle, when it reaches an age 
of 30 Myr. 

For thermal feedback to be effective, the pressure gradient es¬ 
tablished by the heating must be able to perform work on the gas 
(via sound waves, or shocks for supersonic flows) on a timescale 
that is shorter than that required to erase the gradient via radiative 
cooling. By comparing the sound crossing and cooling timescales 


® [Durier & Dalla Vecchia|j?012) show that implementations of kinetic and 
thermal feedback converge to similar solutions in the limit that the cool¬ 
ing time is long. We adopt a thermal implementation here for consistency 
with our AGN feedback implementation, and because numerical tests indi¬ 
cate that it is less susceptible to problems stemming from poor numerical 
sampling of outflows. 


for heated resolution elements, |Dalla Vecchia & Schaye| ( [2012^ de¬ 
rived an estimate for the maximum gas density, nH,t„, at which 
their stochastic heating scheme can be efficient (their equation 18), 




10 cm 


T 


K 


3/2 


106 M© 


- 1/2 


(4) 


where T > ATsf is the post-heating temperature. For simplic¬ 
ity, the cooling is assumed to be Bremsstrahlung-dominated, so the 
value of riH.tc will be over-estimated in the temperature regime 
where (metal-)line cooling dominates (T ^ 10*^ K). S15 noted that 
some stars do form in the EAGLE simulations from gas with den¬ 
sity greater than the critical value, in which case the spurious (nu¬ 
merical) radiative losses are significant. In the case of such losses 
being significant, the energy budget used to reproduce the GSMF 
in the simulation will likely be an overestimate of that required in 
Nature. 

Equation indicates that numerical losses associated with 
stars forming from high-density gas can be mitigated by appeal¬ 
ing to a higher heating temperature, ATsf. However, this is not an 
ideal solution. Eor a fixed quantity of feedback energy per stellar 
mass formed (i.e. constant /th), the probability that a star particle 
triggers a heating event is inversely proportional to ATsf. Based 
on the energy budget described above, the expectation value of the 
number of resolution elements (in this case, SPH particles) heated 
by a star particle is specified by |Dalla Vecchia & Schaye| j2012| 
equation 8) as: 


(Aheat) 


1.3/th 


/ ATsf 
VlO^-s K 


-1 


(5) 


In the regime of ATsf 10*^ ® K, the probability of heating a 
single resolution element becomes small, leading to poor sampling 
of the feedback cycle. We therefore only consider models adopting 
ATsf = 10^ ® K. 

If ATsf is sufficiently high to ensure that numerical losses 
are small, the physical efficiency of feedback can be controlled by 
adjusting /th. This mechanism is a means of modelling the sub¬ 
grid radiative losses that are not addressed by our simple treatment 
of the ISM. Because these losses should depend on the physical 
conditions in the ISM, there is physical motivation to specify /th 
as a function of the local properties of the gas. Primarily, it is the 
freedom to adjust /th that enables the simulations to be calibrated. 


2.5 Black holes and AGN feedback 


Eeedback associated with the growth of BHs is an essential ingredi¬ 
ent of the EAGLE simulations. Besides regulating the growth of the 
BHs, AGN feedback quenches star formation in massive galaxies 
and shapes the gas profiles of their host haloes. The implementa¬ 
tion adopted here consists of two elements, namely i) a prescription 
for seeding galaxies with BHs and for following their growth via 
mergers and gas accretion, and ii) a prescription for coupling the 
radiated energy, liberated by BH growth, to the ISM. The method 
for the former is based on that introduced by |Springel et al.lpOOS i 
and modified by |Booth & Schay^j2009^ and Rosas-Guevara et ak] 
( |2013t , while the method for the latter is similar to that described 
by Booth & Schaye|l 2009^. 


Eollowing Springel et al.| ( [2005] l, seed BHs are placed at the 
centre of every halo more massive than 10*® Mq//i that does not 
already contain a BH. Candidate haloes are identified by running a 
EoE algorithm with linking length b = 0.2 on the dark matter dis¬ 
tribution. When a seed is required, the highest density gas particle 
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is converted into a collisionless BH particle, inheriting the parti¬ 
cle mass and acquiring a subgrid BH mass ttibh = 10® Mq/Zi 
which, for the intermediate resolution simulations considered here, 
is smaller than the initial gas particle mass by a factor of 12.3. Cal¬ 
culations of BH properties are therefore functions of ttibh, whilst 
gravitational interactions are computed using the particle mass. 
When the subgrid BH mass exceeds the particle mass, BH particles 
stochastically accrete neighbouring gas particles such that particle 
and subgrid BH masses grow in concert. BHs are prevented from 
“wandering” out of their parent haloes by forcing those with mass 
< lOOrrig to migrate towards the position of the minimum of the 
gravitational potential in their halo. 

BHs are merged if separated by a distance that is smaller than 
both the BH kernel size, /ibh, and three gravitational softening 
lengths, and if their relative velocity is smaller than the circular ve¬ 
locity at the distance /ibh, v^ei < where /ibh and 

mBH are, respectively, the SPH kernel size and the subgrid mass 
of the most massive BH in the pair. The relative velocity thresh¬ 
old prevents BHs from merging during the initial stages of galaxy 
mergers. 

2.5.1 Gas accretion onto black holes 


that couples to the ISM. In common with the efficiency of feed¬ 
back associated with star formation, /th, the value of et must be 
chosen by calibrating to observations. Because of self-regulation, 
the value of ef only affects the masses of BHs ( [Booth & Schaye| 
|2009| l, which vary inversely with et, and it has little effect on the 
stellar mass of galaxies (provided its value is non-zero). The pa¬ 
rameter et can be calibrated by ensuring the normalisation of the 
observed relation between BH mass and stellar mass is reproduced 
at z = 0. Although implemented as a single heating mode, |Mc-| 
jCarthy et al.| ( |20ll) demonstrate that this scheme mimics quiescent 
‘radio’-like and vigorous ‘quasar’-like AGN modes when the BH 
accretion rate is a small (<C 1) or large (~ 1) fraction of the Ed¬ 
dington rate, respectively. 

S15 demonstrated that the efficiency adopted by OWLS, et = 
0.15, remains a suitable choice at the higher resolution of EAGLE. 
Therefore, a fraction CfCr = 0.015 of the accreted rest mass en¬ 
ergy is coupled to the local ISM as feedback. Each BH maintains 
a “reservoir” of feedback energy, Ebh. After each time step At, 
an energy eterm^ccrf? At is added to the reservoir. Once a BH has 
stored sufficient energy to heat at least one fluid element of mass 
mg, it becomes eligible to heat, stochastically, its SPH neighbours 
by increasing their temperature by ATagn- For each neighbour the 
heating probability is 


The gas accretion rate, rhaccr, is specified by the minimum of the 
Eddington rate. 


47rGmBHmp 
niBdd = “—“— 

CrCrxC 


( 6 ) 


and 


maccr — ruin ^mBondi[(Cs/L0) /Gvisc], mBondi) , (2) 

where muondi is the Bondi-Hoyle ( |1944| l rate applicable to spheri¬ 
cally symmetric accretion. 


47rG^mBHp 

= (C2 + ^2)3/2 ■ 


( 8 ) 


Here nip is the proton mass, ax the Thomson cross section, c the 
speed of light, Cr the radiative efficiency of the accretion disc, and v 
the relative velocity of the BH and the gas. Finally, V,j, is the circu¬ 
lation speed of the gas around the BH computed using equation 16 
of jRosas-Guevara et al.| l [20T3t and Gvisc is a free parameter related 
to the viscosity of a notional subgrid accretion disc. The growth of 
the BH is specified by 


mBH — (1 er)maccr. 


(9) 


We assume a radiative efficiency of e^ — 0.1. The factor 

(Ca/V 0 )®/Gviac by which the Bondi rate is multiplied in equation 
l^is equivalent to the ratio of the Bondi and viscous time scales 
(see jRosas-Guevara et al.|2013 i. The critical ratio of V,f,/cs above 
which angular momentum is assumed to suppress the accretion rate 
scales as . Larger values of Gvisc therefore correspond to a 

lower subgrid kinetic viscosity, and so act to delay the growth of 
BHs by gas accretion and, by extension, the onset of quenching by 
AGN feedback. 


P = 


Ebb 

AeAGNWngb (mg) 


( 10 ) 


where Aeagn is the change in internal energy per unit mass cor¬ 
responding to the temperature increment, ATagn (the parameter 
ATagn is converted into Aeagn assuming a fully ionised gas of 
primordial composition), Angb is the number of gas neighbours of 
the BH and (mg) is their mean mass. The value of Ebb is then 
reduced by the expectation value of the injected energy. The time 
step of the BHs is limited to aim for probabilities P < 0.3. 

Larger values of ATagn yield more energetic feedback 
events, generally resulting in reduced radiative losses (as per equa¬ 
tion]^. However, larger values also make the feedback more inter¬ 
mittent. In general, the ambient density of gas local to the central 
BH of galaxies is greater than that of star-forming gas distributed 
throughout their discs, so a higher heating temperature is required 
to minimise numerical losse^ The EAGLE reference model pre¬ 
sented by S15 adopts ATagn = 10®'® K. In that study an alterna¬ 
tive intermediate-resolution model was also presented (model “AG- 
NdT9”) which, by appealing to ATagn = 10® K, was found to 
more accurately reproduce the observed gas fractions and X-ray 
luminosities of galaxy groups. This higher temperature increment 
was also found to be necessary in high-resolution simulations, since 
they resolve higher ambient densities close to BHs and hence ex¬ 
hibit higher cooling rates. 

The values of the relevant subgrid parameters adopted by all 
simulations featured in this study are listed in Table 


3 CALIBRATED SIMULATIONS 

As discussed by S15 (see their § 2), if subgrid models for energy 
feedback offer an incomplete description of the processes they are 


2.5.2 AGN feedback 

A single mode of AGN feedback is adopted, whereby energy is 
injected thermally and stochastically, in a manner analogous to en¬ 
ergy feedback from star formation (see § |2.4^ . The energy injection 
rate is efCrrhaccrC®, where et is the fraction of the radiated energy 


® BHs are in principle able to inject feedback energy at all times, unlike 
star particles which inject prompt feedback only once. The AGN feedback 
cycle can therefore remain well sampled for higher heating temperatures 
than is the case for star formation feedback, as long as the interval between 
heating events is shorter than a Salpeter time [Booth & Schaye|2009|. 
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Identifier 

Side length 
[ cMpc] 

N 

'Yeos 

[cm 

3] 

/th -scaling 

tfthjmax 

/th,min 

^H,0 

[cm“®] 

nn 

C^visc/^TT 

^kTAGN 
logio [K] 

Calibrated models 

FBconst 

50 

752 

4/3 

Eq. 

3 



1.0 

1.0 



10 ^ 

8.5 

FBo- 

50 

752 

4/3 

Eq. 

3 


2 

'’'dm 

3.0 

0.3 

- 

- 

10 ^ 

8.5 

FBZ 

50 

752 

4/3 

Eq. 

3 


Z 

3.0 

0.3 

- 

- 

10 ^ 

8.5 

Ref (FBZp) 

50 

752 

4/3 

Eq. 

3 


Z,p 

3.0 

0.3 

0.67 

2 / In 10 

10 “ 

8.5 

Ref variations 
eosl 

25 

376 

1 

Eq. 

3 


Z,p 

3.0 

0.3 

0.67 

2 / In 10 

10 “ 

8.5 

eos5/3 

25 

376 

5/3 

Eq. 

3 


Z,p 

3.0 

0.3 

0.67 

2 / In 10 

10 “ 

8.5 

FixedSfThresh 

25 

376 

4/3 

0.1 



Z,p 

3.0 

0.3 

0.67 

2 / In 10 

10 “ 

8.5 

WeakFB 

25 

376 

4/3 

Eq. 

3 


Z,p 

1.5 

0.15 

0.67 

2 / In 10 

10 “ 

8.5 

StrongFB 

25 

376 

4/3 

Eq. 

3 


Z,p 

6.0 

0.6 

0.67 

2 / In 10 

10 “ 

8.5 

ViscLo 

50 

752 

4/3 

Eq. 

3 


Z,p 

3.0 

0.3 

0.67 

2 / In 10 

10 ^ 

8.5 

ViscHi 

50 

752 

4/3 

Eq. 

3 


Z,p 

3.0 

0.3 

0.67 

2 / In 10 

10-2 

8.5 

AGNdTS 

50 

752 

4/3 

Eq. 

3 


Z,p 

3.0 

0.3 

0.67 

2 / In 10 

10 “ 

8.0 

AGNdT9 

50 

752 

4/3 

Eq. 

3 


Z,p 

3.0 

0.3 

0.67 

2 / In 10 

10 “ 

9.0 


Table 1. Parameters that are varied in the simulations. Columns are: the side length of the volume (L) and the particle number per species (i.e. gas, DM) per 
dimension (N), the power law slope of the polytropic equation of state (7eos), the star formation density threshold (n|j), the scaling variable of the efficiency 
of star formation feedback (/th). the asymptotic maximum and minimum values of /tn, the Ref model’s density term denominator (hh.o) and exponent 
(rin) from equation |14| the subgrid accretion disc viscosity parameter (Cvisc) from equation |7] and the temperature increment of stochastic AGN heating 
(ATagn)- The upper section comprises the four models that have been calibrated to reproduce the z = 0.1 GSMF, and the lower section comprises models 
featuring a single-parameter variations of Ref (varied parameter highlighted in bold). All models also adopt riz = 2/ In 10 with the exceptions of FBcr, for 
which the parameter nz is replaced by nx with the same numerical value (see equationfT^, and FBconst, for which the parameter is inapplicable. 


designed to model, are subject to numerical losses, or if the out¬ 
comes of the prescriptions are sensitive to resolution, then the true 
efficiencies of feedback processes cannot be predicted from first 
principles. It was therefore argued that cosmological hydrodynami- 
cal simulations are presently unable to yield ab initio predictions of 
the stellar mass of galaxies, nor the mass of their central BH. Sub¬ 
grid parameters should therefore be calibrated such that simulations 
reproduce desired diagnostic quantities, stellar and BH masses be¬ 
ing germane examples. 

The optimal approach to calibrating subgrid models is not un¬ 
ambiguous, since there can be multiple measurable outcomes that 
are sensitive to the adjustment of subgrid parameters, some or all 
of which might reasonably be considered valid constraints. For ex¬ 
ample, in the case of feedback efficiencies, one might calibrate the 
model to reproduce the velocity and/or mass loading of outflowing 
gas. However, these quantities remain ill-characterised observation- 
ally, and are sensitive to the physical scale on which they are mea¬ 
sured (which is generally not even well known). Reproducing the 
properties of outflows on a particular spatial scale offers no guar¬ 
antee that they are reproduced on other scales, since, for example, 
the interaction of outflows with the circumgalactic medium may be 
inadequately modelled. The choice of calibration diagnostic(s) is 
therefore somewhat arbitrary, but some choices can be more read¬ 
ily motivated. Clearly, it is necessary that any diagnostic be well 
characterised observationally on the scales resolved by the simu¬ 
lation. Perhaps the most elegant example is the star formation law 
which, on the ~ 10^ pc scales we follow in the ISM, can be accu¬ 
rately represented by the Kennicutt-Schmidt relation; as described 
in § |2.2| the free parameters of the subgrid star formation law are 
unambiguously prescribed by observations. In addition, it is desir¬ 
able to confront calibrated models with complementary observa¬ 
tional constraints, to minimise the risk of overlooking modelling 
degeneracies. 

We chose to calibrate the feedback simulations to the z = 0.1 


GSMF, a practice commonly adopted by the semi-analytic galaxy 
formation modelling community. Low-redshift galaxy surveys en¬ 
able the GSMF to be characterised in the local Universe across five 
decades in mass scale, with a precision that is, assuming a univer¬ 
sal IMF, limited primarily by systematic uncertainties in the stellar 
evolution models used to infer the masses jConroy et al.|2009{|Pforr| 
|et al.|201^ but see |Taylor et al.|20n) , peculiar velocity corrections 
for faint galaxies (e.g. Baldry et al.||201 2), and the method used 
to subtract the sky background from bright galaxies (e.g.|Bemardi| 
|et al.|20f3]|Kravtsov et al.|2014| l. An additional motivation for ap¬ 
pealing to the GSMF as the calibration diagnostic is that its re¬ 
production by the simulations is a prerequisite for examination of 
many observable scaling relations. 


Whilst calibrating the simulations, attention was also paid to 
the sizes of galaxies. As with the calibration diagnostic, this choice 
is somewhat arbitrary, but it is readily motivated since the formation 
of unrealistically compact or extended galaxies would limit the util¬ 
ity of the simulations and likely indicate physical and/or numerical 
inaccuracies. The formation of disc galaxies with realistic sizes in 
cosmological simulations has proven to be a non-trivial challenge, 
leading to the identification and rectification of many shortcomings 


of numerical techniques (e.g.|Sommer-Larsen et al. 1999||Ritchie & 

Thomas|20011 Marri & White|20031|0kamoto et al.|20031|Agertz 

et al.||2007 Springel|2010 

Hopkins|2013 120^. In terms of the 

physics of galaxy assembly. 

Navarro & Whit^ 1994J identified the 


transfer of angular momentum from cold, dense clumps of baryons 
to the non-dissipative dark matter residing at the outskirts of galaxy 
haloes as the major concern. In a comprehensive simulation com¬ 


parison project, |Scannapieco et al.| ( [2012^ highlighted that the angu¬ 
lar momentum problem (and also the overcooling problem) remains 
without a consistent solution, and therefore that galaxy sizes cannot 
yet be uniquely predicted, even when the assembly history of their 
parent halo is fully specified. It cannot be assumed that the sizes 
of galaxy discs will be accurately reproduced by simulations, even 
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if they successfully reproduce a suitable calibration diagnostic. For 
this reason, we require a model to reproduce both the GSMF and 
the observed size-mass relation of disc galaxieQat low-redshift, in 
order for the calibration process to be deemed successful. 

The subgrid model governing feedback associated with star 
formation in EAGLE is primarily dependent upon the IMF, /th and 
ATsf- a universal IMF is adopted through out, a nd the adoption of 
a fixed ATsf = 10^'® K was motivated in § 2.4 Therefore, energy 


feedback associated with star formation is calibrated exclusively 
by varying /th, the fraction of the total available energy from type 
II SNe that couples to the ISM. The main effect on the GSMF of 
increasing (decreasing) /th is to lower (raise) its normalisation in 
terms of the comoving number density of galaxies with stellar mass 
below the “knee” of the [SchechterU 1976) function (demonstrated 
in § [4A31 i. 

The subgrid model governing AGN feedback in EAGLE is pri¬ 
marily dependent upon et , Gvisc and ATagn- The retention of the 
AGN feedback efficiency adopted by OWLS (cf = 0.15) was mo¬ 
tivated in § |2.5.2| In contrast to ATsf, whose value is fixed by the 
need to suppress numerical radiative losses and to adequately sam¬ 
ple the feedback process, the freedom to adjust ATagn, which de¬ 
termines the energetics and intermittency of AGN feedback events, 
can be motivated. We therefore explore simulations with different 
AGN heating temperatures but, in terms of the calibration of galaxy 
properties, ATagn only (weakly) affects the stellar mass of galax¬ 
ies with M* > 10 ®^ Mq (it does, however, impact markedly upon 
the properties of the intragroup and intracluster media, e.g. |Le Brun| 
|et al.|2014)|Schaye et al.|20T5) . Eor the purposes of reproducing the 
present-day GSME, AGN feedback is primarily calibrated by vary¬ 
ing Gvisc, which broadly governs the mass scale at which AGN 
feedback becomes efficient. Since this mass scale is weakly depen¬ 
dent upon Gvisc, models can adopt values of this parameter that 
differ by orders of magnitude. 

The calibration of phenomenological components of galaxy 
formation models is a practice that was established during the de¬ 


model^ 

(|Kauffmann et al.|l993 Cole et al. 

1994 

Somerville & 

IPrimack 

1999 Cole et al.|2000 Hatton et al.| 

2003). Semi-analytic 


models are built upon the framework of dark matter halo merger 
trees, so the parameters they adopt for governing feedback must be 
coupled to simplified models for the structure of galaxies and their 
interstellar and circumgalactic media. In contrast, hydrodynamical 
simulations enable feedback properties to be specified, if so de¬ 
sired, based on the physical conditions local to newly formed star 
particles. 

Thus far the capability to specify feedback parameters in this 
fashion has only been exploited by a limited number of groups, 
each using a very similar feedback implementation: kinetically- 
driven winds that are launched outside of the ISM (by decoupling 
the wind particles from hydrodynamic forces), and whose proper¬ 
ties are imposed by specifying their initial velocity and mass load¬ 
ing factor {rj = Mwind/Af*) as a function of the properties of the 
dark matter environment, for example the gravitational potential 
or velocity dispersion. Simulations adopting this implementation 


^ The size evolution of both late- and early-type galaxies will be explored 
in detail in a forthcoming paper (Furlong et al. in prep). 

* Modem semi-analytic models often adopt distribution sampling tech¬ 


niques to efficiently calibrate their many free parametersJKampakoglou 


et al. 

2008 

Henriques et al. 

2009 

Bower et al.|2010||Lu et al.|2012||Mutch 

et al. 

2013 

Henriques et al. 

2013 

2oT4r 


have been used to investigate the establishment of the GSMF jOp- 


penheimer et al.|[2010[ [Puchwein & Springel||201^ |Vogelsberger 


et al.|2013) and the formatio n of disc galaxies similar to the Milky 

Way I Marinacci et al.|2014|, to examine the observed properties of 
some intergalactic absorption systems (jOppenheimer & Dave|2006| 
|2009[fWgelsberge r et al.|2 014b|> and to re produce the Local Group 
satellite population jOkamoto et al.|2010) . 

This implementation is well-motivated since, by temporarily 
decoupling winds and specifying their properties as a function of 
dark matter properties, it affords simulations the best opportunity 
to achieve numerical convergence. However, it precludes the full 
exploitation of the hydrodynamics calculation. The physical prop¬ 
erties of outflows are almost certainly dependent upon the local 
(baryonic) conditions of the ISM, and these properties are avail¬ 
able to use as inputs to subgrid feedback models. Since the phi¬ 
losophy adopted for the EAGLE project is to calibrate the feedback 
scheme, the convergence demands placed upon the adopted subgrid 
models are relaxed, presenting the appealing opportunity to couple 
the value of subgrid parameters (e.g. /th) to the baryonic properties 
of the local ISlvfl 


3.1 Calibrating the star formation feedback efficiency 

The role of /th in shaping the z — 0.1 GSMF is investigated in 
this section. Examination of the previous generation of simulations 
upon which the EAGLE project is based indicates that AGN feed¬ 
back is a necessary ingredient for regulating the growth of mas¬ 
sive galaxies ( |Crain et'aL||2009[ [Schaye et al.||201(){ |Haas et ^ 
|2013i^ and establishin g the gas-phase properties of galaxy groups 
i McC^hy et al.|20T0l|2011[|Le Brun et al.|2014) . Four calibrated 
simulations are explored, each featuring energy feedback associ¬ 
ated with both star formation and the growth of BHs. All four sim¬ 
ulations adopt ATagn = 10®'® K, and each features a constant 
value of Gvisc that is allowed to differ between simulations such 
that, when combined with the function specifying /th, the result¬ 
ing z = 0.1 GSMF features a break close to the observed scale of 


In three of the models, asymptotic efficiencies of /th.max = 3 
and /th.min = 0.3 are used. Values greater than unity can be mo¬ 
tivated physically, since there are sources of energy feedback other 
than type II SNe, and indeed such sources are often invoked in sim¬ 
ulations of galaxy formation, for example stellar winds and radia¬ 
tion pressure ([Stinson et al.|2013{|Hopkins et al.| 2014[ l or cosmic 
rays jJubelgas et al.||2008[ [Booth et al.||20'T3 1. However the pri¬ 
mary motivation for appealing to /th > 1 here is to offset numer¬ 
ical losses that result from the finite resolution of the simulations. 
There are two means by which finite resolution introduces artificial 
losses. The first, as discussed in §2.4| (equation[4l(, being that there 
exists a maximum density above which stochastic thermal feedback 
is inefficient, because the pressure gradient established by feedback 
is erased by radiative cooling before it is able to exert mechanical 
work on the gas. The second stems from the inability of large cos¬ 
mological simulations to model the formation of the earliest gen¬ 
eration of stars. As discussed by |Haas et al.| ( [2013a) , this means 
that the first generation of galaxies that form in simulations do so 


® SIS introduced the nomenclature “weak convergence” to describe the 
consistency of simulation outcomes in the case that subgrid parameters are 
recalibrated when the resolution is changed, as opposed to “strong conver¬ 
gence” in the case of holding the parameters fixed. 
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'o9io Tqm [K] 

3 4 5 6 7 8 



'°9io Z [Z©] 


Figure 1. The fraction of the energy budget due to type II SNe feedback 
that is used for thermal heating, /th, in the FBconst, FBct, FBZ and Ref 
models. The FBconst model, represented by the dashed grey line, adopts 
/tjj = 1, independent of local conditions, /tn declines smoothly as a func¬ 
tion of Tdm in the FBcr model (equation |12) , and as a function of Z in the 
FBZ model (equation 1 1 3j . The upper axis is aligned and scaled such that 
both FB(t (upper axis) and FBZ (lower axis) are described by the dark blue 
curve (no physical coirespondence between Tom nnd Z is implied by this 
alignment). The Ref model adds a density dependence to FBZ (equation 
|14| , such that for stars forming from gas with nj] < riH.o ths /th function 
is shifted to lower values (e.g. cyan curve for nn = nH,o/3) and vice versa 
(e.g. red curve for njj = Srij] o)- 


within haloes that have not been subject to feedback, and hence ex¬ 
hibit unrealistically high gas fractions and star formation efficien¬ 
cies. Appealing to /th > 1 for these galaxies partially compensates 
for this unavoidable shortcoming. 


FBconst 

The simplest model injects into the ISM a fixed quantity of energy 
per unit stellar mass formed, independent of local conditions. This 
value corresponds to the total energy liberated by type II SNe, i.e. 
/th = 1- The adopted subgrid viscosity parameter for BH accretion 
is Cvisc = 27r X 10®. Although the injected energy is independent 
of local conditions, scale-dependent macroscopic radiative losses 
can nonetheless develop self-consistently, for example due to dif¬ 
ferences in the metallicity (and hence cooling rate) of outflowing 
gas, the ram pressure at the disc-CGM interface, or the depth of 
the potential well. This model therefore provides a baseline against 
which it is possible to assess the degree to which the overall phys¬ 
ical losses need to be established by calibrating losses on subgrid 
scales. Because /th = 1 represents the uncalibrated case, there 
is no reason to expect that FBconst will reproduce the observed 
z = 0.1 GSMF. However, we will see later that it does do so, but 
fails to reproduce the observed sizes of disc galaxies. 


FBa 


Okamoto et al. 2010[ Oppenheimer et al. 2010 Puchwein & 


Springel|2013[|Vogelsberger et al.pof^ Khandai et al.|2014[ Vo 

gelsberger et al. |2014at . However, because these studies all adopt 

kinetically-driven outflows imposed outside of the ISM, and spec¬ 
ify the properties of the outflows by imposing a mass loading and 
initial wind velocity, our implementation is significantly differ¬ 
ent. Rather than specifying the properties of the outflows, sim¬ 
ilar behaviour on macroscopic scales is sought as an outcome 
of the stochastic heating implementation, by calibrating the effi¬ 
ciency, /th, as a function of ctom- The latter is the square of the 
3-dimensional velocity dispersion of dark matter particles within 
the smoothing kernel of a star particle at the instant it is born, and 
is a proxy for the characteristic virial scale of the star particle’s 
environment. 


^ ^ 5 / ^DM \ 2 

3fc ^ ^^VlOOkms-1/ 


( 11 ) 


For simplicity, we assume at all times the mean molecular weight 
of a fully ionized gas with primordial composition, jj, = 0.591. 
The fit to the z = 0.1 GSMF for this model is achieved with a 
slightly higher subgrid viscosity for BH accretion than is the case 
for FBconst, Cvisc = 271 X 10®. 

Since the properties of star formation-driven outflows are 
linked to the state of local dark matter only via gravitational forces, 
no physical motivation for /th(TbM) is sought. The adopted func¬ 
tional form simply maximises the feedback efficiency (by minimis¬ 
ing putative subgrid radiative losses) in low-mass galaxies, whilst 
reducing the feedback efficiency in more massive counterparts, 
where the conversion of gas into stars is known to be most ef¬ 
ficient (e.g. |Eke et al.||2005| |Behroozi et al.||2013| |Moster et al.| 
|2013| >. At higher masses still, AGN feedback is assumed to domi¬ 
nate the regulation of star formation, so a low star formation feed¬ 
back efficiency for high-dispersion environments is reasonable. The 
adopted functional form of /th is a logistic (sigmoid) function of 
logio Tdm, 


/th — /th,inin T 


/th ,max - /th ,n 


1 + 


( ^DM \ ' 
10 ® K ) 


( 12 ) 


shown in Figure[2 (dark blue curve, corresponding to the upper 
x-axis). The function asymptotes to /th.max and /th.min in the 
limits Tdm <C 10® K and Tdm 10® K, respectively, and 
varies smoothly between these limits about Tdm = 10® K (or 
(Tdm — 65 km s“®). The parameter nx > 0 controls how rapidly 
/th varies as the dark matter “temperature” scale deviates from 
10® K. The rather unnatural value nx = 2/ In 10 ~ 0.87 follows 
from an early implementation of the functional form adopted in the 
feedback routine; an exponent of unity would yield similar results. 


FBZ 

Adjusting the subgrid radiative losses with the metallicity of the 
ISM assigns a physical basis to the functional form of /th- Phys¬ 
ical losses associated with star formation feedbaclp^ are likely to 
be more significant when the metallicity is sufficient for cooling 
from metal lines to dominate over the contribution from H and He. 
For temperatures 10® K < T < 10^ K, characteristic of outflow¬ 
ing gas in the simulations, the transition is expected to occur at 


This model adopts the popular convention of prescribing feed¬ 
back parameters according to local conditions inferred from neigh¬ 
bouring dark matter particles ^Oppenheimer & Dave|20()6l |2008[ 


A metallicity dependence for tf (§ |2.5.2f is not explored, because metals 
are not expected to dominate the radiative losses at the higher temperatures 
associated with AGN feedback. 
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Z ~ 0.1 Zq jWiersma et al.||2C)09b^ . This qualitative behaviour 
is captured by the same functional form as equation [T^ replacing 
(TbM,nT,10® K) with (Z,nz,0.1 Zq) to obtain, 


fth — /th.min 


fth .max — fth .n 


1 + 


[O.lZgJ 


(13) 


where Zq = 0.0127 is the solar metallicity ( [Allende Prieto et"^ 
20011 and nz = ut = 2/In 10. This function corresponds to 
the dark blue curve and the lower x-axis in Figure[T] with fth 
asymptoting to /th.max and /th.min in the limits Z < 0.1 Z© and 
Z 3> 0.1 Z©, respectively. Since galaxies tend to follow a tight re¬ 
lation between their mass and metallicity, the feedback efficiency 
is, as in FBct, greatest for low-mass galaxies. Moreover, because 
metallicity characteristically decreases with redshift at fixed stel¬ 
lar mass, the feedback efficiency is weighted towards early cos¬ 
mic epochs. This helps to partially decouple the growth of galaxies 
from the growth of their parent halo, which is thought to be a nec¬ 
essary condition for reproducing the observed number density evo¬ 


lution of low-mass galaxies in a CDM cosmogony (e.g. Weinmann 


et al. 2012[|Flenriques et al.|2013]|Mitchell et al.|20T4{ [White et al. 


2014 1 . The subgrid viscosity parameter for AGN is the same as per 
FBcr, Cvisc = 2n X 10^. 


Ref (equivalently, FBZp) 


A significant fraction of the star particles in the FBcr and FBZ mod¬ 
els form at densities greater than nn.tc. the critical density above 
which feedback energy is quickly radiated away (equation|^. The 
feedback associated with these high-density star formation events 
is therefore numerically inefficient. The consequences of this over¬ 
cooling are explored later in § |4.1[ The overestimated losses can be 
compensated by introducing a density dependence in the expres¬ 
sion for fth- 

fth ,max — fth ,n 


fth — /th,min T 


( Z birth \ 

V "H.O ) 


(14) 


where nn, birth is the density of a gas particle at the instant it is 
converted into a star particle. The feedback efficiency therefore in¬ 
creases with density at fixed metallicity, whilst respecting the orig¬ 
inal asymptotic values. The choice of nn.o = 0.67 cm“® was 
guided by a suite of small test simulations, which also indicated 
that the adoption of rin = nz is sufficient to reproduce the z = 0.1 
GSMF. The effect of the additional density term is illustrated in 
Figure where it can be seen that for nn = nn.o the functional 
form adopted by Ref is identical to that of FBZ, but the fth curve 
is shifted to lower (higher) values for stars forming from lower 
(higher) density gas. Only the shift to higher fth for higher density 
gas (at fixed metallicity) is required to offset numerical losses, but 
for simplicity we choose to include the shift to lower efficiency at 
low density that is implied by the adopted function. Such a density 
dependence may also have a physical basis: because the star for¬ 
mation law has a supra-linear dependence on surface density, the 
feedback energy injection rate per unit volume increases with den¬ 
sity. At fixed density, a higher energy injection rate corresponds to 
higher temperatures and longer cooling times. Therefore, feedback 
associated with clustered star formation is expected to lead to lower 


radiative losses (e.g.|Heiles|1990| Creasey et al. 

2013 [Krause et al.| 

|2013||Nath & Shcheldnov|2013 Roy etal.|2013 

Keller etal.|2014J, 


and vice-versa. The density-dependent feedback efficiency adopted 
here (equation |14|> ensures that, when integrating over all star for¬ 
mation events in the simulation, the mean and median values of fth 


remain close to unity (they are 1.06 and 0.70, respectively, for Ref- 
L0100N1504). The subgrid viscosity parameter for BFl accretion is 
lower in Ref than in the other calibrated models, Gvisc = 2?! x 10°. 


4 RESULTS 

This section begins with an examination of the calibrated simula¬ 
tions. In §4.2| simulations featuring single parameter variations to 
the reference model are briefly introduced, and the impact of the 
changes on the resulting galaxy population are explored. 


4.1 Examination of the calibrated simulations 

We begin by examining the models calibrated to reproduce the 
2 : = 0.1 GSMF. Besides the GSMF, the evolution of the comoving 
stellar mass density of each simulation is examined, as are the sizes 
and specific star formation rates (SSFRs) of galaxies alz = 0.1. As 
discussed by S15, the GSMF and galaxy sizes were used for the cal¬ 
ibration process, so are not presented as predictions. Although the 
model parameters were not adjusted in order to improve the corre¬ 
spondence between the simulated and observed SSFRs, the former 
were inspected throughout the calibration process, and hence the 
predictions of the SSFRs are not “blind”. The aim of this exercise 
is to examine how different implementations of physical processes 
impact upon the properties of galaxies and their environments. As 
part of this procedure, the conditions of the ISM from which stars 
are born are also explored. 


4.1.1 The galaxy stellar mass function 


The 2 ; = 0.1 GSMFs of the four calibrated simulations, FBconst 
(red curve), FBct (green), FBZ (cyan) and Ref (dark blue), run us¬ 
ing L050N0752 initial conditions, are shown in the left-hand panel 
of Figure]^ The right-hand panel also shows the GSMF of the ref¬ 
erence model at intermediate resolution in a L = 100 cMpc vol¬ 
ume (Ref-L100N1504, red curve), which was introduced by S15, 
and in a smaller L = 25 cMpc realisation (Ref-L025N376, green), 
which is used in §4.2| as a baseline against which to compare sev¬ 
eral single parameter variations of the reference model. Maintain¬ 
ing the convention established hy SI5, curves are drawn with dot¬ 
ted lines where galaxies are less massive than 100 (initial mass, 
trig) baryonic particles, as resolution tests (presented in S15) indi¬ 
cate that sampling errors due to finite resolution are significant in 
this regime. At the high-mass end, curves are drawn with dashed 
lines where the GSMF is sampled by fewer than 10 galaxies per 
(A login M* = 0-2) bin. Data points represent the GSMF in¬ 
ferred from the Sloan Digital Sky Survey ( |Li & White|2009 SDSS, 
filled circles) and the Galaxy and Mass Assembly ( [Baldry et al.| 
|2012| GAMA, open circles) survey. In both cases, the data have 
been adjusted for consistency with the value of the Hubble param¬ 
eter adopted by the simulations. GAMA measurements on scales 
M* < 10® M© are drawn as lower limits, since the data are ex¬ 
pected to be incomplete at the typical surface brightness associated 
with these galaxies ( [Baldry et al.|2012l ). 

In the range of stellar masses where the mass resolution and 
volume of the simulation enable a robust measurement of the 
GSMF (10® < M* < 10^^ M©), the number density of galaxies at 
fixed stellar mass produced by each of the four calibrated models 
is consistent with the observational data to < 0.31 dex. This pre¬ 
cision is comparable to the systematic uncertainty associated with 
spectrophotometric techniques for inferring galaxy stellar masses. 
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Figure 2. Left: The 2 = 0.1 GSMF of the calibrated L050N0752 simulations, FBconst (red), FBcr (green), FBZ (cyan) and Ref (dark blue). Curves are drawn 
with dotted lines where galaxies are comprised of fewer than 100 star particles, and dashed lines where the GSMF is sampled by fewer than 10 galaxies 
per bin. Data points show measurements with Icr error bars from the SDSS l |Li & White|2009l filled circles) and GAMA i jBaldry et al.|2Ql^ open circles) 
surveys. The simulations each reproduce the observed number density of galaxies at fixed mass to < 0.31 dex, a precision unprecedented for hydrodynamical 
simulations. Right: To illustrate convergence as the simulation volume is varied, the Ref model at intermediate resolution is shown in volumes of L = 25, 50 
and lOOcMpc. The Ref-L050N0752 (dark blue) and Ref-LI0ON1504 (red) GSMFs are consistent for M* < lO^^'^ Mq. The GSMF of Ref-L025N0376 is 
consistent with that of its larger counterparts for M* < 10^’^ Mq, but samples large-scale structures poorly owing to its small volume, imprinting “wiggles” 
on the GSMF. 


indicating that a more precise reproduction of the GSMF may be 


that a more precise rep: 
T3(e .g.|Conroy et al.p 


2009 


Pforr et al. 


2012 I. This de- 


unwarrantei 

gree of consistency with observational data is typical of that associ¬ 
ated with semi-analytic galaxy formation models. The reproduction 
of 2: = 0.1 GSMF by multiple cosmological hydrodynamical sim¬ 
ulations featuring feedback efficiencies governed by such distinct 
schemes is unprecedented in the literature (see also Fig. 5 of S15). 

The detailed confrontation of the EAGLE reference model 
with observational data presented by S15 was based on the Ref- 
L100N1504 simulation. Since running multiple L100N1504 sim¬ 
ulations is, at present, computationally prohibitive, the calibrated 
models have been run with L050N0752 initial conditions. It is 
therefore necessary to confirm that the Ref-L050N0752 simu¬ 
lation yields a GSMF that is consistent with that of its L = 
lOOcMpc counterpart. Comparison of the Ref-L050N0752 and 
Ref-L100N1504 curves in the right-hand panel of Figure indi¬ 
cates that the GSMFs of these simulations are consistent to a high 
precision (< 0.063 dex at fixed stellar mass) over the range for 
which both simulations are well sampled (M* < 10^^'® Mq). 
An L = 50cMpc volume is therefore sufficient to capture the 
effects of subgrid physics on all but the most massive galax¬ 
ies seen in observational surveys of the local Universe. The Ref- 
L025N0376 simulation tracks its larger counterparts on scales of 
M* < 10 ®'® Mq, but lacks the volume required to sample more 
massive scales precisely, as is clear from the “wiggles” imprinted 
on the GSMF. 

We do not explore resolution convergence here; S15 demon¬ 
strated the strong convergence behaviour of the 2 = 0.1 GSMF 
for the reference model, and established that recalibration of the 
subgrid parameters enables competitive and well-understood weak 


The GSMF is also impacted upon by other systematic effects, such as 
completeness and extinction corrections, and background subtraction. 


convergence behaviour for a broad range of observational diagnos¬ 
tics. 


4.1.2 Galaxy sizes 

Following |McCarthy et al.H2012^ , we characterise the morphology 
of galaxies by fitting Sersic profiles to their projected, azimuthally- 
averaged surface density profiles. The size of a galaxy is then 
equated to its effective radius, Rso, the radius enclosing 50 per¬ 
cent of the stellar mass when the profile is integrated to infinity. 
The scaling of this quantity with stellar mass for disc galaxies is 
shown in Figure]^ As in |Shen et al.| ( |200^ , whose size measure¬ 
ments from SDSS data are overplotted, disc galaxies are defined to 
be those with Sersic indices Ua < 2.5. The binned median is plotted 
for each calibrated simulation. Dashed lines are used where the me¬ 
dian is sampled by fewer than 10 galaxies, and a dotted linestyle is 
used for bins corresponding to stellar masses < 600mg; this scale 
was shown by S15 to be the minimum for which size measurements 
are robust. In the regime between the limits of sufficient resolution 
and adequate galaxy sampling, the la (i.e. 16**' to 84*^*' percentile) 
scatter about the median of Ref is shown as a blue shaded region. 

The reference model tracks the observed relation closely, with 
the medians of the simulation and the observational measurements 
being offset by 0.1 — 0.2 dex. This offset is comparable to the 
systematic offset between the median of those data and the me¬ 
dian sizes measured by |Baldry et al.|j2012^ based on GAMA ob¬ 
servations of blue galaxies (defined using an r-band magnitude- 
dependent colour threshold). Also shown in Figure are the size 
measurements presented by S15 for the larger Ref-LIOONI 504 vol¬ 
ume (yellow curve) demonstrating that galaxy sizes are unaffected 
by the volume of the simulation, as expected. In contrast to Ref, 
the FBconst, FBct and FBZ models are inconsistent with the ob¬ 
served size-mass relation. Galaxies with M* > 10^® Mq cease to 
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Figure 4. Face-on (left) and edge-on (right) projections, at z = 0.1, of the galaxy formed within the same M 200 ~ 10^^ Mq halo in the four simulations that 
were calibrated to reproduce the present-day GSMF. Each image subtends a field of view of 100 pkpc, and is a composite of SDSS u-, g- and r-band emission 
maps generated with the radiative transfer softwai'e SKIRT. The overall size of the optical envelope of the galaxy is similar in each case, but the distribution and 
dynamics of the stars differ markedly between Ref and the other models. In the latter, the galaxy forms an unrealistically compact bulge component at early 
epochs and exhibits too little ongoing star formation (blue-coloured concentrations) in the extended disc. Consequently, only in Ref does the galaxy exhibit an 
effective radius, R 50 , that is consistent with the observed size-mass relation for disc galaxies. 


follow the observed relation between size and mass, and become 
much too compact, with median sizes only a few times the gravi¬ 
tational softening scale (cprop = 0.7 pkpc). For this reason, these 
models are not considered satisfactory when applying the EAGLE 
calibration criteria. We note that the measurements of both IShenI 
|et al.|j200^ and |Baldry et al.|p012^ are based on r-band photom¬ 
etry, and that isophotal radii are sensitive to the band in which they 
are measured. [Lange et al.|p015) l recently exploited the multiband 
photometry of the GAMA survey to assess the magnitude of this 
sensitivity, concluding that the observed size of disc galaxies of 
mass M* = 10 ^° Mq is typically only ~ 10 percent smaller in 
the Lfg-band (the best proxy for the true stellar mass distribution) 
than in the r-band. We can therefore be confident that our con¬ 
clusions here are unlikely to be strongly affected by, for example, 
attenuation by interstellar dust. Detailed tests of mock observations 
derived by coupling EAGLE to the SKIRT radiative transfer algo¬ 
rithm ( jBaes et al.|2011) will be presented in a forthcoming study 
(Trayford et al. in prep). 

The calibrated simulations adopt identical initial conditions, 
enabling galaxies to be compared individually as well as statisti¬ 
cally. In Figure|^ the galaxy that forms within the same dark matter 
halo is shown for each of the calibrated simulations at z = 0.1. The 
halo was selected at random from those with mass in the Ref simu¬ 
lation 8 X 10^^ M 0 < M 200 < 2 X 10^^ Mq. In each simulation, it 
therefore hosts a galaxy whose stellar mass (M* ~ 2 x 10 ^° M©) 
corresponds to the minimum of the size-mass relation exhibited by 
the FBconst, FBcr and FBZ simulations. Each image subtends a 
field of view 100 pkpc on a side, and is a composite comprised 
of monochromatic SDSS u-, g-, and r-band emission maps. The 
maps were generated with SKIRT (jBaes et ^[2011|), which con¬ 


siders the photometric properties of the stellar populations and the 
estimated dust distribution, the latter being inferred from the pre¬ 
dicted metallicity of the ISM. The same mapping between physical 
flux and pixel luminosity is adopted in each panel. The galaxy is 
shown face-on (left-hand panels) and edge-on (right-hand panels), 
oriented about the angular momentum vector of the star particles 
comprising the galaxy. 

Visual inspection indicates that the outer envelope of the 
galaxy, corresponding to an r-band surface brightness of ~ 
28 mag arcsec”^, is similar in each simulation. However, the dis¬ 
tribution of stars within that radius differs markedly between Ref 
and the other models, and this strongly influences the effective ra¬ 
dius. In the FBconst, FBct and FBZ simulations, the galaxy forms 
a massive, compact bulge that dominates the overall stellar dis¬ 
tribution. In the Ref simulation the star-forming disc component, 
seen clearly in the face-on images as blue-coloured concentrations 
distributed over all radii, comprises a greater fraction of the mass. 
Based on a dynamical decomposition similar to the orbital circu¬ 
larity method of jAbadi et al.| ( [2003[ (, the bulge-to-total ratio of the 
galaxy in the FBconst, FBcr, FBZ and Ref simulations is 0.47, 0.43, 
0.50 and 0.30, respectively. 


Many studies in the literature conclude that the ability of a 
hydrodynamical simulation to reproduce, approximately, the low- 
redshift GSMF (or, in the case of ‘zoom’ simulations, the M* — 
Mhaio relation as inferred by e.g. subhalo abundance matching) 
to be a metric of success, without comparing to the observed size 
of galaxies (e.g. Okamoto^^e^^lJj201BjJOppenheirner^e£^alJ2010l 


Munshi et al.|2013[|Puchwein & Springel|2013[|Vogelsberger et al.j 


2013 1 . However the reproduction here of the observed 2 ; = 0.1 
GSMF with a number of models that yield unrealistically compact 
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Figure 3. The sizes, at z = 0.1, of disc galaxies in the four simulations that 
were calibrated to reproduce the present-day GSMF. Sizes from the Ref- 
L100N1504 (presented by S15) are also shown (yellow curve), to demon¬ 
strate consistency between volumes. Size, i? 5 o, is defined as the half-mass 
radius (in proper kpc) of a Sersic profile hi to the projected, azimuthally- 
averaged stellar surface density profile of a galaxy, and those with Sersic 
index ris < 2.5 are considered disc galaxies. Curves show the binned me¬ 
dian sizes, and are drawn with dotted lines below a mass scale of 600mg, 
and a dashed linestyle where sampled by fewer than 10 galaxies per bin. 
The 1(7 scatter about the median of Ref is denoted by the blue shaded re¬ 
gion. The solid and dotted grey lines show the median and la scatter of 
sizes for Ug < 2.5 galaxies inferred from SDSS data by |Shen et al.H2003) , 
whilst grey data points and error bars show sizes of blue galaxies infen'ed 
by |Baldry et al.|j2012) from GAMA data. Only the Ref model successfully 
reproduces the observed GSMF and galaxy sizes at z = 0.1. 

galaxies, highlights the importance of simultaneously calibrating 
models with observational diagnostics that are complementary to 
the GSMF. In Sections |4.1.3| |4.1.4| and |4.1.5| we turn to the ex¬ 
amination of diagnostics that were not considered during the cal¬ 
ibration process. It is shown that models that yield unrealistically 
compact galaxies also fail to reproduce the observed star formation 
history and present-day specific star formation rate of the galaxy 
population. We demonstrate that the formation of compact galax¬ 
ies is a consequence of numerical radiative losses becoming severe 
in high density gas, thus artificially suppressing the efficiency of 
energy feedback. 

4.J.3 Comoving stellar mass density 

By construction, the calibrated simulations yield similar volumetric 
stellar mass densities at z = 0.1. The evolution of stellar mass den¬ 
sity in the four simulations can differ, however, because the history 
of energy injection from feedback varies between the models. Fig¬ 
ure]^ shows the evolution of the comoving, instantaneou stellar 
mass density in the four calibrated L050N0752 simulations. The 
total comoving density of stars in each simulation is shown; in a 
companion paper, [Furlong et ar]p014| l excluded diffuse intraclus¬ 
ter light (by considering only those stars within 30 pkpc of galac¬ 
tic centres) to mimick observational measurements, and recovered 

Stellar evolution mass loss by star particles is accounted for. 


Figure 5. The evolution of the comoving stellar mass density of the cali¬ 
brated EAGLE models. Data points correspond to measurements from sev¬ 
eral observational surveys spanning 0 < z < 4 (see text for details). Al¬ 
though the models each broadly reproduce the observed GSMF at z = 0.1, 
they exhibit markedly different star formation histories, and the FBconst 
model in particular forms too much stellar mass at early times. Taking the 
observations at face value, only the Ref model is broadly consistent with 
the observational constraints for z < 2. 


densities in the Ref-L100N1504 simulation that were lower by ap¬ 
proximately 20 percent for z < 1. 

Data points represent the comoving stellar mass density in¬ 
ferred from a number of complementary observational analyses. 
Where necessary, the data have been adjusted to adopt the same 
IMF and Hubble parameter as the simulations. The filled and open 
circles at z ~ 0.1 represent the integration of the SDSS JLi &| 
jWhitejgOO^ and GAMA jBaldry et al.][20T2) GSMFs shown in 
Figure [^ respectively. Diamonds represent measurements over the 
redshift interval 0.1 < z < 0.9 inferred from a combined sam¬ 
ple of SDSS and PRIMUS data presented by [Moustakas et al.j 
l |20I3^ , triangles represent measurements over the redshift inter¬ 
val 0.2 < z < 4 inferred from Ultra VISTA data by jMuzzin et al.j 
(|20I3f, and squares represent measurements over the redshift inter¬ 
val 0.625 < z < 2.25 inferred from ZFOURGE data by |Tomczak| 
jet al.j p014^ . Data from surveys that overlap in redshift interval 
are shown in order to illustrate, broadly, the degree of systematic 
uncertainty and field-to-field variance in the measurements. 

The stellar mass densities of the four simulations are consis¬ 
tent to < 0.1 dex at z = 0.1. The FBconst simulation, however, 
forms stars too rapidly at early epochs. The stellar mass density 
inferred from observations at z ~ 1 is in place in this simula¬ 
tion prior to z = 2, whilst the evolution at intermediate redshifts 
(1 < z < 3) is, by necessity, then too weak. The models that al¬ 
low the star formation feedback efficiency to vary as a function of 
the local environment track the observed build up of stellar mass 
more accurately, since they typically inject more energy (per unit 
stellar mass formed) into star-forming regions in low-mass galax¬ 
ies (which dominate at high redshift), than is the case for FBconst. 
However, the FB(t and FBZ models remain inconsistent with the 
observational measurements at z > 1, and only the Ref model 
broadly reproduces these out to z ~ 2. 
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4.1.4 Specific star formation rates 


The markedly different evolution of the comoving stellar mass den¬ 
sity in the calibrated simulations is indicative of similar differ¬ 
ences between the models in terms of the star formation rates of 
star-forming galaxies z = 0.1. The relation between the SSFR 
and stellar mass at z = 0.1 is adopted as the diagnostic 
with which to compare the models; the evolution of the SSFR of 
Ref is explored by [Furlong et al.H2014^ . Figure shows the me¬ 
dian SSFR of star-forming galaxies, as a function of stellar mass, 
at 2 = 0.1 for the four calibrated simulations. As per Figure]^ 
the median curve is drawn with a dashed linestyle when sampled 
by fewer than 10 galaxies per bin. The diagonal dotted line indi¬ 
cates the SSFR corresponding to 10 star-forming gas particles at a 
gas density of nn = 0.1 cm“®; at SSFRs below this limit, sam¬ 
pling limitations are significant and median curves are drawn with 
a dotted linestyle. In the regime between the limits of sufficient 
resolution and adequate galaxy sampling, the la scatter about the 
median is shown for Ref as a blue shaded region. The horizon¬ 
tal dashed line denotes the threshold SSFR of 10“^ Gyr“^, that 
was adopted to separate star-forming from passive galaxies. Data 
points represent SSFR measurements of star-forming galaxies in¬ 
ferred from SDSS- Stripe 82 data ([Gilbank et al.|2010[ squares) and 
the GAMA survey (|Bauer et al.|2013 circles). 


S15 demonstrated that, at M* ~ 10^^ Mq, the Ref- 

L100N1504 simulation exhibits SSFRs very similar to those ob¬ 
served, but at lower masses the simulated SSFRs are systematically 
lower than observed by up to ~ 0.3 dex. This remains the case for 
Ref-L050N0752. However, as shown by S15, the discrepancy for 
low mass galaxies is much smaller for the Recal-L025N0752 sim¬ 
ulation, indicating that our high-resolution simulations are better 
able to reproduce the star formation properties of low-mass galax¬ 
ies. The FBZ model behaves similarly to Ref. The FBconst and 
FB(t models also exhibit SSFRs similar to those observed in mas¬ 
sive galaxies, but at 10 ® M© they are ~ 0.7dex lower than ob¬ 
served. 


In the framework of equilibrium galaxy formation models (e.g 
Finlator & Dave|2008[ [Sdiaye et al.|2010[[bave et al.|2012[rMi-| 

tra et al.||2014| ), the gas inflow rate onto galaxies is balanced by 

the combined sinks of star formation and ejective feedback, and 
the specific inflow rate (at fixed redshift) which is a weak func- 


tion of halo mass 

Dekel et al.|20091 

Fakhouri et al.|20101|van de| 

[Voort et al.[[2011| 

Correa et al.|[20H 

J. Since the calibrated sim- 


ulations each broadly reproduce the observed GSMF, galaxies of 
fixed stellar mass occupy similarly massive haloes in each case: at 
2 = 0 . 1 , the median halo mass associated with galaxies of stellar 
mass Mi, = 10® M© in the Ref simulation is offset from that of the 
FBconst simulation by < 0.1 dex. This leaves differences in the 
mass reaccretion rate and the efficiency of preventive feedback as 
prime candidates for establishing an offset in the present-day SSFR 
of low-mass galaxies. 

It is indeed likely that the reaccretion of ejected gas is sensi¬ 
tive to the details of the feedback (e.g. |Oppenheimer et al.|2010| 
[Brook et al.|2014| > and almost certainly plays a role in shaping the 
present-day SSFR of galaxies, particularly so at the mass scale cor¬ 
responding to M*. We intend to explore this process in detail in 
a forthcoming study (Crain et al. in prep). The efficiency of pre¬ 
ventive feedback is, by construction, a distinguishing feature of 
the four calibrated models, and one that is simple to explore. Ex¬ 
amination of the median value of /th associated with star forma¬ 
tion events over the gigayear preceding 2 = 0.1 for galaxies of 
Mi, ~ 10® Mq shows marked differences: for the FBconst, FBcr, 



log,o M. [Mq] 

Figure 6. The median specific star formation rate, M*/M*, of star-forming 
galaxies as a function of stellar mass at 2 = 0.1, of the calibrated EAGLE 
simulations. The diagonal dotted line corresponds to the SSFR of 10 gas 
particles at rijj = 0.1 cm“^, below which sampling effects are significant 
and medians are drawn with a dotted linestyle. The Icr scatter about the 
median of Ref is denoted by the blue shaded region. The dashed horizontal 
line denotes the SSER separating star-forming and passive galaxies. Data 
points with error bars correspond to the median and Icr scatter of the SSER 
of star-forming galaxies inferred from SDSS-Stripe 82 data by jGilbanket al. [ 
j2010[ squares) and GAMA data by[Bauer et al.[j2013[ circles). 


FBZ and Ref models, the values are 1, 1.04, 0.58 and 0.35, re¬ 
spectively. The star formation rate required to produce sufficiently 
strong outflows from star formation feedback scales inversely with 
these efficiencies, and thus the Ref model correspondingly exhibits 
the highest SSFR at 2 = 0.1. 


4.1.5 The birth conditions of stars 

The properties of simulated galaxies are clearly sensitive to the 
adopted functional form of /th. Galaxy sizes, which encode infor¬ 
mation related to the state of the gas from which stars were bom, 
are the clearest discriminator of the models explored here, indicat¬ 
ing a connection between the stucture of the ISM and the efficacy 
of feedback. 

When star particles are bom, we record the density of their 
parent gas particle, enabling an examination of the physical condi¬ 
tions of the gas from which all stars in the simulations were bom. 
The EAGLE simulations treat star-forming gas as a single-phase 
fluid, therefore the SPH density of star-forming particles can be 
considered as the mass-weighted average of the densities of cold, 
dense molecular clouds and of the warm, ionised medium with 
which they maintain a pressure equilibrium. Pressure is therefore 
a more physically meaningful property of star-forming gas in the 
simulations, and it is possible to recover the birth pressure of stars 
from their birth density under the reasonable assumption that their 
parent gas particle resided on the Jeans-limiting pressure floor at 
the time of conversiorF^ 

Figure]^ shows the differential distribution of ISM pressures 


Particles within 0.5 dex of the temperature associated with the pressure 
floor (see §[2.2| are eligible to form stars. 
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Figure 7. Differential distribution of the pressure of fluid elements (i.e. 
SPH particles) at the instant they were converted to star particles, for all 
star particles formed prior to 2 = 6 (dotted curves), 2 = 2 (dashed) 
and 2 = 0.1 (solid). The upper axis shows the ratio of the correspond¬ 
ing density to the critical density for numerically efficient feedback (equa¬ 
tion!^. At high redshift, most stars form from metal-poor gas with a pres¬ 
sure corresponding to the maximum star formation threshold (equation j^, 
P/k ~ 10®'® Kcm“®, though in the FBconst model a significant frac¬ 
tion of stars has already formed from gas with nn ^H,tc ■ As the ISM 
becomes enriched, the threshold drops and stars are able to form at lower 
ISM densities. However at 2 < 6 spurious numerical losses in the FBconst, 
FB(t and FBZ simulations result in the rapid build up of a second peak at 
high-pressure. The development of this peak is suppressed in the Ref sim¬ 
ulation owing to the greater efficiency of feedback at higher density (for 
fixed metallicity), which ensures that the majority of stars are able to yield 
numerically efficient feedback. 


(normalised by Boltzmann’s constant, k) at the instant of their for¬ 
mation, of star particles formed in the calibrated simulations prior 
to 2 = 6 (dotted curves), 2 = 2 (dashed curves) and 2 = 0.1 (solid 
curves). The upper axis indicates the corresponding ratio of the gas 
particle density to the critical density for which stochastic thermal 
heating associated with star formation is efficient (nn.tci equation 
1^. Examination of this ratio affords us a means by which to test 
for numerical overcooling on an event-by-event basis. At high red- 
shift, most stars form from low-metallicity gas, and are hence sub¬ 
ject to a high-star formation density threshold (equation j^. The 
maximum of this threshold is nn = 10cm“®, corresponding to 
P/k ~ 10®'® K cm“® for our choice of Jeans limiting equation of 
state. Many stars form from gas with pressures close to this thresh¬ 
old value. 

A significant fraction of stars also form from higher-pressure 
gas in the FBconst simulation prior to 2 = 6. The formation of 
stars from gas with nn > leads to artificial radiative losses. 

As discussed in the fact that the first galaxies whose forma¬ 
tion can be captured by the simulations are associated with haloes 
that have not been subject to feedback, means that they exhibit arti¬ 
ficially high gas fractions and star formation efficiencies. This ini¬ 
tial problem has the potential to set in train a cycle of overcool¬ 
ing: the artificially rapid initial formation of stars over-enriches the 
ISM with efficient coolants, promoting further cooling losses and 
enabling dissipation to higher densities. Stars subsequently form¬ 
ing from this gas yield numerically inefficient thermal feedback 



Figure 8. The birth pressure of star particles, as a function of galactocen- 
tric radius, for galaxies with stellar mass 10 < logj^Q M* / Mq < 10.5 
at 2 = 0.1. Solid curves denote the binned median, and the blue shaded 
region corresponds to the Icr scatter about the median of Ref. The dotted 
vertical line shows the gravitational softening scale. The highest pressures 
correspond to stars formed in galactic centres, so numerically inefficient 
feedback in the FBconst, FBcr and FBZ models leads to the formation of un¬ 
realistically massive central spheroids. Conversely, the suppression of high 
birth densities in Ref yields effective radii that are consistent with obser¬ 
vations. The convergence of the birth pressure profiles for r > 4pkpc 
indicates that differences in the heating rate of gas (i.e. /th) between the 
models are confined to galactic centres. 

(because nn > nH.t^), so the gas fraction and star formation ef¬ 
ficiency of the halo remain artificially high. An initial numerical 
shortcoming therefore has the potential to trigger unrealistic phys¬ 
ical losses that themselves promote further numerical losses. This 
cycle can lead to a strong overestimate of the severity of radiative 
losses. 

The adoption of = 3 for stars forming in low-velocity 
dispersion (FBcr) and low-metallicity (FBZ, Ref) environments ef¬ 
fectively eliminates the initial phase of the problem; at 2 = 6, by 
which time the simulations comprise tens of thousands of star parti¬ 
cles, the fraction of stars formed from high-pressure gas is small for 
the FBct, FBZ and Ref simulations. As the ISM becomes enriched 
with metals, the typical star formation threshold drops and a peak in 
the distribution of birth densities develops at P/k ~ 10® '® K cm“® 
(riH ~ 0.3 cm“®) in each of the calibrated simulations. This in¬ 
jection of additional energy into nascent galaxies is, however, in¬ 
sufficient to arrest the onset of subsequent numerical losses. At 
2 < 6, the birth pressure distribution of the FBct and FBZ sim¬ 
ulations (in addition to that of FBconst) develops a second peak 
at P/k ~ 10^ ®Kcm“®, corresponding to nn ~ 250 cm“®. At 
such high density, resolution elements heated to 10^'® K cool be¬ 
fore they can expand, rendering thermal feedback numerically in¬ 
efficient. The fraction of stars formed from gas with nu > tiH.tc 
offers a simple estimate of the severity of this numerical overcool¬ 
ing; for stars formed prior to 2 = 0.1 the fractions for the FBconst, 
FBct, FBZ and Ref simulations are 0.55, 0.56, 0.61 and 0.25, re¬ 
spectively. 

The suppression of the high-pressure peak in the distribution 
of birth pressures exhibited by the Ref simulation is achieved by 
adding a density-dependence to the star formation feedback effi- 
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ciency function (equation 1 14^. It ensures that star formation feed¬ 
back remains efficient when stars form from relatively dense gas, 
preventing the build-up of the highest pressures within the ISM. A 
density dependence of this sort can also be motivated on physical 


grounds (e.g.|Heiles|1990| 

Creasey et aL|2013||Krause et al.|2013| 

|Nath & Shchekinov|2013[ 

Roy et al.|20131|Keller et al.|2014J, but 


our main motivation is to combat the numerical problems described 
above. 


Figure shows the median birth pressure of stars as a func¬ 
tion of their galactocentric radius at 2 = 0 . 1 , for galaxies with stel¬ 
lar mass in the interval 10^° < M*/M 0 < This interval 

corresponds to the mass scale at which the difference in sizes be¬ 
tween Ref and the FBconst, FBcr and FBZ models is greatest. The 
blue shaded region shows the Icr scatter about the median of Ref, 
and the vertical dotted line denotes the gravitational softening scale 
of 0.7 pkpc. The correpondence between numerically inefficient 
feedback and the formation of unrealistically compact galaxies is 
clear: the highest birth pressures are exhibited by stars residing in 
the centres of galaxies. The reduction of star formation from gas 
with P/k > 10®'® Kcm“® in the Ref simulation therefore prefer¬ 
entially suppresses the formation of compact galactic bulges, and 
enables the formation of galaxies with effective radii that are con¬ 
sistent with the observed size-mass relation. The median birth pres¬ 
sure profiles converge above a characteristic radius of r ~ 4 pkpc, 
indicating that the difference in heating rates between the simula¬ 
tion (due to differences in the adopted functional form of fth) is 
only significant within the centres of galaxies. 



Figure 9. The ratio of the stellar mass to halo mass of central galaxies, as 
a function of stellar mass and normalised by the cosmic baryon fraction, in 
L0025N0376 simulations adopting 7eos = 1 (eosl, green), 4/3 (Ref, dark 
blue) and 5/3 (eos5/3, red). Curves show the binned median ratios, and are 
drawn with dotted lines below the mass scale corresponding to 100 baryonic 
particles, and a dashed line where sampled by fewer than 10 galaxies per 
bin. The Icr scatter about the median of Ref is denoted hy the blue shaded 
region. Dark and light grey lines represent the abundance matching relations 
of |Behroozi et al.|f2013) and |Moster et al.|p013^ , respectively. 


4.2 Variations of the Reference model 

The Ref model demonstrates that it is possible to calibrate a model 
that satisfactorily reproduces the GSMF and the observed sizes of 
galaxies at 2 = 0.1. It is important to quantify the sensitivity of 
the outcomes of this model to variation of its key subgrid parame¬ 
ter 0 This is achieved by conducting a series of simulations within 
which the value of a single parameter is varied from that adopted 
by Ref, as listed in the lower section of Table [T] Parameters gov¬ 
erning the ISM, star formation and the efficiency of star forma¬ 
tion feedback are tested using relatively inexpensive L025N0376 
simulations. Those governing the AGN feedback are tested with 
L050N0752 simulations, since the effects of changing these pa¬ 
rameters are most clearly imprinted upon the properties of massive 
galaxies and their environments. The effects of reasonable changes 
can vary markedly from parameter to parameter, so we focus here 
only on properties of the galaxy population that shift significantly 
from those of the corresponding Ref simulation. Simulations where 
the variation does have a significant effect are likely to yield a 
galaxy population that is no longer an accurate representation of 
the observed Universe. 

4.2.1 ISM equation of state 

The simulations adopting isothermal ( 7 eos = 1, “eosl”) and adia¬ 
batic ( 7 eos = 5/3, “eos5/3”) equations of state are examined in this 
section. Steeper slopes, i.e. higher values of 7 eos, yield higher ISM 
pressures as the volume density is increased, and hence a shorter 
gas consumption time at fixed density (assuming a star formation 
law with n > 1 , which is the case here) and a suppression of the 

The effects of changes to the hydrodynamics and time-stepping schemes 
will be explored in a forthcoming companion paper (Schaller et al. in prep). 


development of high ISM densities. For 4/3 < 7 eos < 2, the Jeans 
conditions indicate that gas is able to collapse without fragmenting. 
Steeper slopes inhibit collapse, whilst shallower slopes may lead to 
artificial fragmentation (e.g. |Bate & Burkert|1997||Schaye & Dalla| 
|Vecchia|2008) . 

The effects of varying 7 eos were previously explored in ide¬ 
alised disc galaxy simulations by |Schaye & Dalla Vecchia 1 2008 1 
and in the OWLS cosmological simulations by |Schaye et al.' 20101 
an d |Haas et al.| ( |2013b^ . These studies concluded that, as long as 
the index was maintained between values of unity and 5/3, the pri¬ 
mary effect of varying 7 eos is upon the visual appearance of the 
disc, with the higher pressures associated with stiffer equations of 
state yielding smoother gas distributions and larger scale heights. 
Otherwise, differences due to variations in 7 eos were found to be 
minimal, and this was attributed to the fact that the pressure-law 
implementation of star formation adopted by OWLS (and EAGLE, 
§ | 2 . 2 ^ reproduces the observed star formation law, irrespective of 

'7’eos- 

Figure|^shows the ratio of the stellar mass of central galaxies 
to the mass of their host halo, normalised by the cosmic average 
baryon fraction (Ob/Om = 0.157), and plotted as a function of 
stellar mass. The dark and light grey lines represent the median stel¬ 
lar mass-to-halo mass ratios inferred from the abundance match¬ 
ing algorithms of |Behroozi et al.| ( [2013| l and |Moster et al.|f2013| l, 
respectively. For low-mass galaxies, the conclusion that the prop¬ 
erties of galaxies are largely insensitive to 7 eos is consistent with 
EAGLE. Changing 7 eos does, however, have a significant impact 
on the most massive galaxies since, in contrast to the OWLS ref¬ 
erence model, the EAGLE reference model includes treatments of 
the growth of BHs and the associated AGN feedback. In the regime 
that the local sound speed is large compared to the relative velocity 
of the central BH and the surrounding gas (Cs S> v), the accre- 
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tion rate (equation onto the BH is inversely proportional to the 
cube of the ISM sound speed. For a polytropic equation of state, 
the sound speed Cg oc P/p oc so it is clear that a stiffer 

equation of state for the ISM will yield a greater sound speed at 
fixed density, and so reduce the accretion rate onto the BH. 

Adopting an isothermal (adiabatic) equation of state leads to 
galaxies with M* > 2 x 10^° Mq exhibiting a median BH mass 
that is approximately 0.3 dex larger (smaller) than that of Ref. This 
difference in BH mass, and the associated difference in the rates of 
BH accretion and AGN feedback, translates into significant differ¬ 
ences in the mass of gas converted to stars by z = 0.1 in halos 
of mass M 200 s/ 10^^ M©, as shown in Figure]^ The efficiency 
of AGN feedback is therefore sensitive to the assumed polytrope 
7 eos, and variation of the latter would require a recalibration of 
the subgrid BH viscosity parameter, Gvisc, to recover the observed 

GSMI0 

4.2.2 Star formation threshold 

A simulation adopting a constant density threshold for star for¬ 
mation of riH = 0.1 cm“® (“FixedSfThresh”), as used by de¬ 
fault in the OWLS and GIMIC projects, is used to examine the 
role of the dependence of the star formation density threshold, 
n^{Z) on metallicity, as implemented in Ref. The use of constant 
and metallicity-dependent density thresholds for star formation was 
also investigated using the OWLS simulations by |Scha ye et al.| 
( |2010^ and |Haas et al.| ( [2013b^ . As in that study, the low-redshift 
galaxy population exhibits no significant differences between sim¬ 
ulations run with a constant density threshold of njj = 0.1 cm“® 
and the metallicity-dependent threshold described by equation 
therefore we do not show comparisons of these simulations here. 
This can be understood by appealing to self-regulation arguments 
( [Schaye et al.|2010) . 


4.2.3 Star formation feedback efficiency 

There is a growing consensus, based on analyses of cosmological 
simulations, that galaxy formation and evolution is governed pri¬ 
marily by the regulation of the supply of gas to the ISM, rather than 


the detailed physics of interstellar gas itself (e.g.|Rasera & Teyssier| 


2006| Bouche et al.|2010 Schaye et al.|2010[ 

Dutton et al.|2010[ 

Dave et al.||2012 Haas et al.||2013a|b Vogels 

berger et al.||20131(. 


In general, it is feedback that regulates this fuel supply, so it is 
reasonable to expect that changing the efficiency of star formation 
feedback will have a significant impact on many characteristics of 
the galaxy population. This can be gauged by inspection of Figure 
|10| which shows the effect, on a broad range of galaxy scaling rela¬ 
tions, of scaling the /th function adopted by the Ref model (equa- 
tion[^ by factors of 0.5 (WeakFB) and 2 (StrongFB). The asymp¬ 
totic efficiencies of these models are therefore = (0.15, 0.6) 
and = (1-5, 6.0), respectively. The median efficiency associ¬ 
ated with all star particles formed prior to 2 = 0.1 in the WeakFB, 
Ref and StrongFB L0025N0376 simulations are 0.38, 0.63 and 
1 .22, respectively. 

As in the previous plots, binned medians are shown (WeakFB: 
green curve, Ref: dark blue curve, StrongFB: red curve), drawn 
with a dotted linestyle below the estimated resolution limit, specific 


In principle, it is possible to construct a model of BH growth that, like 
the star formation implementation (equation|^, is independent of 7eos ■ This 
would require that the accretion rate be specified as a function of pressure. 


to each diagnostic, and a dashed linestyle where there are fewer 
than 10 galaxies per bin. Where appropriate, the Ict scatter about 
the median of Ref is shown as a blue shaded region. The panels 
are labelled (a) to (i) in row-by-row and top-to-bottom order, show¬ 
ing (source of observational data in parentheses): (a) the GSMF 
( [Baldry et al.|2012) ; (b) the stellar mass to halo mass ratio of cen- 
tral galaxies jBehroozi et al.|2013J; (c) the stellar half-mass radius 
( |Shen et al.|2003> , (d) the maximum circular velocity, which we 
adopt as a proxy for the Tully-Fisher relation ( |Avila-Reese et alT] 
|2008|l, (e) the sp ecific star formation rate of star-forming galax¬ 
ies (Bauer et al.||2013|l; (f) the fraction of galaxies that are pas- 
( Moustakas et al.|2073l l; (g) the central BH mass dMcCom^ 


& Ma||2013^; (h) the metallicity of star-forming gas dZahid et ~ 


2014 1 , and (i) the metallicity of stars ([Gallazzi et al.||2005 1 . For 


more details concerning the observations and their comparison with 
the simulations, see S15. 

Inspection of the panels highlights that relatively small 
changes (i.e. factor 2) to the value of /th adopted by Ref has a 
dramatic impact upon the properties of the 2 = 0.1 galaxy popula¬ 
tion. Panel (a) shows that a lower (higher) star formation feedback 
efficiency corresponds to a greater (lesser) abundance of galaxies 
with masses below M*, the scale corresponding to the “knee” of 
the |Schechter| ( [T97^ function. The cause of this shift of number 
densities is clear from inspection of the stellar mass to halo mass 
relation, (b). The form of the median relation is similar in all three 
cases, but with a significant normalisation offset: in the WeakFB 
(StrongFB) model, galaxies acquire a lower (higher) characteristic 
value of /th throughout their growth. In the framework of equilib¬ 
rium models, this requires a higher (lower) SFR to produce suffi¬ 
ciently strong outflows to balance the inflow of gas and, by 2 = 0.1, 
leads to galaxies of a fixed stellar mass becoming associated with 
less (more) massive dark matter haloes, with respect to Ref. Since 
the dark matter halo mass function is particularly steep (e.g. |Jenk-| 
|ins et al.|2001| >, even a small change to the relation between stellar 
mass and halo mass impacts significantly upon the number density 
of galaxies at a fixed stellar mass. Even on more massive scales, 
the adoption of very efficient star formation feedback dramatically 
reduces the abundance of galaxies, although on these scales the 
adoption of inefficient star formation feedback does not result in a 
commensurate increase in the abundance of galaxies. As explored 
below, in the absence of efficient star formation feedback, BHs sim¬ 
ply grow more massive (at fixed halo mass) in order to liberate the 
energy required to achieve self-regulation. 

The shift of the typical halo mass associated with galaxies of 
fixed stellar mass, as the feedback efficiency is varied, impacts sig¬ 
nificantly upon galaxy scaling relations, as we explore below. In 
general, the following panels highlight the importance of populat¬ 
ing haloes with galaxies of the correct stellar mass, whether “by¬ 
hand” in empirical models such as the halo occupation distribution 
or abundance matching methodologies, or as a result of the imple¬ 
mented baryon physics in semi-analytic and hydrodynamic simula¬ 
tions. 

The connection between feedback and the sizes of galaxies 
was highlighted in Sections [4 L2|and [4.L5| Consistent with that 
discussion, the impact of the star formation feedback efficiency 
upon the sizes of galaxies is evident in panel (c). As per Figure 
only galaxies for which the best-fitting Sersic profile has in¬ 
dex ris < 2.5 are considered. At fixed stellar mass, the WeakFB 
(StrongFB) model yields smaller (larger) galaxies, with respect to 
Ref. The adoption of WeakFB leads to significant overcooling (in 
a physical sense), and results in artificially compact galaxies. In 
the regime of less significant overcooling, the efficiency of feed- 
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Figure 10. Galaxy scaling relations at z = 0.1 of L0025N0376 simulations adopting star formation feedback decribed by equation [l4l (Ref, dark blue), and 
this function scaled by factors of 0.5 (WeakFB, green) and 2 (StrongFB, red). With the exception of panels (a) and (f), the curves show binned medians, and 
the blue shaded region corresponds to the la scatter about the median of Ref. Curves are drawn with dotted lines below the estimated resolution limit of each 
diagnostic, and with dashed lines where sampled by fewer than 10 galaxies per bin. The panels show (a) the galaxy stellar mass function; (b) the stellar mass 
to halo mass ratio of central galaxies; (c) the stellar half-mass radius, (d) the maximum circular velocity, (e) the specific star formation rate of star-forming 
galaxies; (f) the fraction of galaxies that are passive; (g) the central BH mass; (h) the metallicity of star-forming gas, and (i) the metallicity of stai's. The source 
of the observational data {dark grey) for each panel is given in §|4.2.3| 
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back still impacts upon galaxy sizes: feedback preferentially ejects 
the lowest angular momentum gas in galaxies (e.g. |Brook et al.| 
|201 1[[20T2^ , so the heating (and concomitant ejection) of more gas 
per unit stellar mass formed in the StrongFB model with respect to 
Ref increases the median angular momentum of the ISM gas that 
remains to form stars, leading to the formation of more extended 
galaxies. 

Panel (d) shows the maximum circular velocity of haloes, 
kmax = max ( \/GM{< r)/r), as a function of the stellar mass of 
their central galaxies. This is a close relative of the|Tully & Fisher| 
\\911) relation exhibited by disc galaxies, so once again only those 
galaxies with < 2.5 are plotted. Changing the efficiency of star 
formation feedback impacts upon both the zero-point and the slope 
of the relation. The zero-point is affected since, as in (h), for a broad 
range of halo masses, the primary effect of changing the efficiency 
is to change the stellar mass that is associated with a given halo 
mass (see also |McCarthy et al.|2012| >. To first order Knax is a rea¬ 
sonable proxy for halo mass, so changing the star formation feed¬ 
back efficiency essentially shifts the stellar mass associated with a 
fixed V)nax, i.e. translates the median curves left-to-right. However, 
the slope also changes because star formation feedback impacts 
most significantly upon low-mass galaxies, as is also clear from 
inspection of (b). Moreover, in high-mass galaxies weak feedback 
can result in an increase of Vmax due to the formation of compact 
bulges. 

The SSFR as a function of stellar mass is shown in panel (e), 
and highlights that the adoption of weaker (stronger) star formation 
feedback leads to lower (higher) SSFR at fixed stellar mass. At first 
glance this may appear to contradict the interpretation of the rela¬ 
tions presented in Figure]^ (see § |4.1.4[ >, where we concluded that 
the weaker feedback at low redshift in Ref (with respect to FBconst, 
for example) enabled low-mass galaxies to achieve a higher SSFR. 
There is, however, an important distinction with respect to the cal¬ 
ibrated simulations: galaxies of a fixed stellar mass in the WeakFB 
and StrongFB models are not associated with haloes of a similar 
mass to their counterparts in Ref. The association of galaxies of 
fixed stellar mass with less (more) massive haloes in the WeakFB 
(StrongFB) model leads to them experiencing lower (higher) infall 
rates (both from the cosmological accretion of intergalactic gas, 
and the reaccretion of ejected material). Increasing the efficiency 
of feedback enables galaxies to balance a fixed net inflow rate at a 
lower SFR, but this is insufficient to compensate for the increased 
inflow rate that stems from being associated with a more massive 
halo. A related symptom of the different SSFRs seen when varying 
the feedback efficiency, is a significant shift of the fraction of galax¬ 
ies that are passive (SSFR < 10“^ Gyr“^) at 2 : = 0.1, as shown 
in panel (f). The impact of changing the star formation feedback 
efficiency is not a simple shift in the normalisation of passive frac¬ 
tion as a function of stellar mass, since approximately half of all 
galaxies in the stellar mass range examined are classified as pas¬ 
sive in WeakFB. In contrast, very few galaxies with stellar mass 
M* ~ 10^° Mq are passive in StrongFB, whilst at higher masses 
the StrongFB model converges on the relation realised by Ref. The 
adoption of weak star formation feedback results in a dramatically 
greater passive fraction for two reasons. Firstly, galaxies consume 
a greater fraction of their low entropy gas for star formation, and do 
so at early times, reducing the reservoir of cold gas available for star 
formation at the present epoch. The second, and more important ef¬ 
fect, is that changing the efficiency of star formation feedback also 
impacts upon the relationship between the masses of galaxies and 
their central BHs, as shown in panel (g). 

The significantly higher (lower) normalisation of the Mbh — 


M* relation in WeakFB (StrongFB) with respect to Ref may ap¬ 
pear, at first glance, to be counter-intuitive. Using OWLS simula¬ 
tions feat uring the|Booth & Scha ye| ( |2009^ implementation of AGN 
feedback, [Booth & Schaye ( 20 10^ concluded that the mass of BHs 
is established primarily by the mass of their host dark matter halo, 
with a secondary dependence upon the concentration of the halo. 
Naively, one might therefore expect that galaxies of fixed stellar 
mass would host less (more) massive BHs, with respect to Ref, 
in WeakFB (StrongFB), since they reside in lower- (higher-) mass 
haloes. However, this reasoning neglects another consequence of 
changing the star formation feedback efficiency: in order to achieve 
self-regulation, BHs in the WeakFB (StrongFB) case must com¬ 
pensate for the lower (higher) star formation feedback efficiency 
by liberating more (less) AGN feedback energy. Since we adopt a 
fixed feedback efficiency for AGN, er = 0.15, BHs can only adjust 
their energy output by growing at different rates and hence yielding 
Mbh — M* relations that are offset from that of Ref at 2 = 0.1. 
Note that in the WeakFB case, no upturn is seen in the Mbh — M* 
relation at intermediate mass scales, since the lack of stellar feed¬ 
back enables BHs hosted by low-mass galaxies to reach the mass at 
which their growth becomes regulated, or even quenched, by AGN 
feedback. 

Finally, panels (h) and (i) show the metallicities of the ISM 
(specifically, star-forming particles) and stars, respectively, as a 
function of stellar mass. At the low- and intermediate-mass scales, 
the simple shift in the normalisation of the relations as a function 
of the adopted stellar feedback efficiency indicates that the outflows 
driven by this mechanism eject metal-rich gas from galaxies, pre¬ 
venting galaxies from acquiring the metallicity expected in the sim¬ 
ple “closed box” chemical evolution scenario. The convergence of 
the median relations for the three models at high metallicity in¬ 
dicates that the ejection of metal-rich gas in massive galaxies is 
less efficient than at low masses, and/or is not driven primarily by 
star formation feedback (see also e.g. E. Rossi et al.|2007[ [D ave| 

|et al.|20lT| [Obreja et al.|2014[ [Creasey et al.|2015^ . The observa- 

tion of a similar anti-correlation between baryonic mass and metal- 
loss in local galaxies is widely perceived as convincing evidence 
for the ubiquity of outflows, and of their efficiency as a mechanism 
to transport metals away from the ISM (e.g. |Tremonti et al.|2004[ 
IZahid et al.|M4t . 

In summary, the efficiency of feedback associated with star 
formation has a strong effect on a broad range of galaxy scaling re¬ 
lations. This is primarily, but not exclusively, because factor of two 
changes to the efficiency significantly alter the relationship between 
stellar mass and halo mass. The changes seen in the properties of 
the 2 = 0.1 EAGLE galaxy population as a result of such changes 
are, for several key scaling relations, at a level that is significantly 
greater than the observational uncertainty. 


4.2.4 Subgrid accretion disc viscosity 

The parameter Cvisc is related to the inverse of the viscosity of a no¬ 
tional subgrid accretion disc, and has two effects. Firstly, it governs 
the angular momentum scale at which the accretion switches from 
the relatively inefficient viscosity-limited regime to the Bondi- 
limited regime (with both cases being subject to the Eddington 
limit). Secondly, it governs the rate at which gas transits through the 
accretion disc when the viscosity-limited regime applies. A higher 
subgrid viscosity, which corresponds to a lower value of the viscos¬ 
ity parameter Gvisc, therefore leads to an earlier onset of the dom¬ 
inance of AGN feedback, and a greater energy injection rate when 
in the viscosity-limited regime. Figure [TT| shows the 2 = 0.1 seal- 
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Figure 11. Galaxy scaling relations at z = 0.1 in L0050N0752 simulations adopting Cvisc = Svr X 10^ (ViscLo, green), Cvisc = 27 r X 10® (Ref, dark 
blue), and Cvisc = Stt X 10“^ (ViscHi, red). Curves show binned medians, and the blue shaded region comesponds to the Icr scatter about the median of Ref. 
Curves are drawn with dotted lines below the estimated resolution limit of each diagnostic, and dotted lines where sampled by fewer than 10 galaxies per bin. 
The panels show (a) the Mbh — JkI^, relation, (b) the stellar mass to halo mass ratio of central galaxies; (c) the effective radius of disc galaxies as a function 
of stellar mass. Details of the observational measurements (dark grey) are provided in §|4.2.3| 


ing relations that are most significantly affected as Cvtac is varied 
by factors of 10^ (ViscLo) and 10“^ (ViscHi), respectively, from 
the value adopted by Ref. These are (a) the Mbh — M* relation, 
(b) the stellar mass to halo mass ratio of cenhal galaxies, and (c) 
the effective radius as a function of stellar mass. 

The Mbh — M* relation exhibits three distinct regimes: an ini¬ 
tial slow and inefficient growth from the seed mass (10® h~^ M©), 
followed by a period of rapid growth of the BH towards the 
observed high-mass scaling relation and, finally, steadier growth 
along the scaling relation. The slow initial growth stems from sev¬ 
eral physical causes. Growth by mergers with seed mass BHs is 
inefficient, simply because the integrated mass of seeds encoun¬ 
tered by any given BH is small (see Figure 4 of [Booth & Schaye| 
|2009| >. Growth by gas accretion is also initially inefficient, firstly 
because the Bondi rate scales strongly with BH mass (oc M|h)> 
and secondly because the rotational support of the ISM in low- 
mass galaxies is sufficient to maintain the angular momentum of 
gas close to the BH above the accretion threshold introduced by 
[Rosas-Guevara et al.H2013) . The significance of the latter is clear 
from inspection of panel (a), which shows that the characteristic 
mass at which BHs begin to grow efficiently is sensitive to Cvisc- 
This scale is M* ~ 10^° Mq in Ref, whilst in the ViscLo and Vis¬ 
cHi models it is shifted to M* ~ 10^° ® M© and M* ~ 10® M©, 
respectively. 

Once gas accretion is efficient, BHs grow rapidly, because the 
feedback liberated by their growth is initially unable to regulate the 
accretion rate. Once sufficiently massive, however, BHs become 
able to regulate, or even quench, their own growth by gas accretion. 
In this regime, the gas that does accrete onto BHs has relatively low 
specific angular momentum ( jRosas-Guevara et al.||2013l >. There¬ 
fore the stellar mass at which BHs arrive on the scaling relation 
is less sensitive to Cvisc than the stellar mass at which BHs begin 
to accrete efficiently. The slope of the Mbh — M* relation in the 
regime of efficient growth is hence sensitive to Cvisc, being signifi¬ 
cantly steeper (shallower) than Ref for ViscLo (ViscHi). However, 
the three simulations do not converge to the same relation, indicat¬ 
ing that accretion onto the BH remains (partially) viscosity-limited 
in even the most massive BHs. In the most massive galaxies, poorly 


sampled by the L = 50 cMpc volumes used here, the growth of 
BHs is likely dominated by mergers. 

The role of the subgrid viscosity as a means to calibrate the 
simulations is clear from panel (b). By shifting the stellar mass 
scale at which BHs begin to self-regulate, the assumed viscosity 
effectively controls the halo mass scale at which AGN feedback 
becomes significant. The amplitude of the baryon conversion effi¬ 
ciency, and the halo mass scale at which the peak occurs, are there¬ 
fore sensitive to the viscosity; appealing to lower (higher) viscosity 
increases (decreases) both the amplitude and halo mass scale of the 
peak stellar mass to halo mass ratio (see also [Rosas-Guevara et al.j 
|20T3l l. 

In the previous section, we concluded that the size of low-to- 
intermediate mass galaxies is sensitive to the efficiency of feedback 
associated with star formation (i.e. /th). Star formation in more 
massive galaxies is regulated primarily by AGN feedback, so their 
sizes are more sensitive to the parameters governing AGN feed¬ 
back, as shown in panel (c). Delaying the onset of efficient AGN 
feedback with a low subgrid viscosity results in the delivery of gas 
to the ISM being countered by star formation feedback alone. A 
greater fraction of the stars in massive galaxies form from particu¬ 
larly dense, metal-rich gas that, as discussed in § |4.L5| is not con¬ 
ducive to the realisation of numerically efficient feedback. Some 
degree of overcooling is therefore to be expected, a symptom of 
which is the formation of a massive, compact bulge component 
that reduces the effective radius of the galaxy. As for adjustments 
to /th, the change in the typical halo mass associated with galaxies 
of a fixed stellar mass also affects the sizes of galaxies, but here the 
effect is weaker, and limited to massive galaxies. 

4.2.5 AGN heating temperature 

The crucial role of the AGN heating temperature in EAGLE was ex¬ 
plored in part by S15, who presented results from a L0050N0752 
simulation adopting ATagn = 10® K. They concluded that the 
higher heating temperature, which yields more energetic but less 
frequent AGN feedback episodes, was necessary to reproduce the 
gas fractions and X-ray luminosities of galaxy groups. Recently, [Lej 
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Figure 12. As per Figure fTT] but varying the AGN heating temperature, ATagn. instead of the subgrid viscosity. The curves correspond to ATagn = 10® K 
(AGNdTS, green), ATagn = 10® ® K (Ref, dark blue), and ATagn = 10® K (AGNdT9, red). 


|Brun et ^ ( [2014^ used the cosmo-OWLS simulation suite to con¬ 
duct a systematic examination of the properties of galaxy groups in 
response to variation of ATagn - They too concluded that a higher 
heating temperature yields more efficient AGN feedback. 

In Figure[^ we again show the Mbh — relation, baryon 
conversion efficiency and the size of galaxies at z — 0.1, this 
time for simulations adopting ATagn = 10® K (“AGNdTS”) 
and ATagn = 10® K (“AGNdT9”), besides Ref which adopts 
ATagn = 10®'® K. The trends appear qualitatively similar to 
those resulting from changes to the subgrid viscosity, but there are 
important differences. Variation of the heating temperature does not 
affect the galaxy stellar mass scale at which BHs begin to grow to¬ 
wards the observed Mbh — M* relation. Therefore, the halo mass 
scale at which the baryon conversion efficiency peaks is also inde¬ 
pendent of the heating temperature. 

However, once a central BH is sufficiently massive to impact 
upon its environment, the macroscopic efficiency (which should not 
be confused with the subgrid efficiency, Cf) of the associated AGN 
feedback is sensitive to the adopted heating temperature. The less 
frequent but more energetic heating events associated with a higher 
heating temperature are more effective at regulating star formation 
in massive galaxies, and as a result the peak baryon conversion ef¬ 
ficiency is higher (lower) in the AGNdTS (AGNdT9) model, with 
repect to Ref. Since the mass of BHs in the simulation is established 
to first order by halo mass jBooth & Schaye|2010^ , the shift in halo 
mass associated with galaxies of stellar mass M* > 10^® Mq es¬ 
tablishes the offset in the Mbh — M* relation above this mass scale. 

The behaviour of the size-mass relation in response to varia¬ 
tion of the AGN heating temperature, panel (c), is consistent with 
that resulting from variation of the star formation feedback effi¬ 
ciency (Figure [Tol> and the subgrid viscosity (Figure pT|. The re¬ 
duced efficacy of AGN feedback when a lower heating temperature 
is adopted leads to the formation of more compact galaxies, since a 
greater fraction of stars are able to form from gas whose density is 
higher than the critical density for numerically efficient feedback, 
n-H.tc (equation|^. 


5 SUMMARY AND DISCUSSION 

This study presents results derived from simulations from the “Evo¬ 
lution and Assembly of GaLaxies and their Environment” (EA¬ 
GLE) project. The simulations adopt cosmological parameters that 
are guided by recent results from the Planck satellite, were realised 
with a force resolution of 0.7 pkpc for z < 3 (better at higher red- 
shift), and adopt a baryonic particle mass of 1.8 x 10® M©. They 
were conducted with a version of the GADGETS software that in¬ 
corporates modified implementations of SPH, time stepping and 
subgrid models governing cooling, star formation, stellar evolu¬ 
tion, feedback from star formation, gas accretion onto and merg¬ 
ers of BHs, and AGN feedback. These simulations complement the 
three introduced in a companion paper by SI5, which also included 
simulations with higher resolution than the intermediate-resolution 
simulations presented here. The EAGLE suite of simulations also 
includes very-high resolution ‘zoom’ simulations of cosmological 
environments similar to the Local Group ( [Sawala et al.|2014a|b^ . 

S15 discussed the implications of the crucial role played by 
subgrid routines in cosmological simulations, concluding that cos¬ 
mological simulations are at present unable to estimate the effi¬ 
ciency of feedback processes from first principles, and thus cannot 
predict stellar and BH masses. The optimal recourse is therefore 
to calibrate the parameters of subgrid routines to ensure that sim¬ 
ulations reproduce well-characterised observables, for example the 
stellar mass function of galaxies. The introduction of several new 
simulations here enables a clear illustration of the motivation for 
the parametrisation adopted by the EAGLE reference model. The 
additional simulations also enable an exploration of the sensitivity 
of the outcomes of the reference model to the variation of its key 
subgrid parameters. 

The EAGLE simulations each adopt the stochastic thermal 
feedback scheme of |Dalla Vecchia & Schaye| ( |2012| ), whereby fluid 
elements are stochastically heated to a pre-specified heating tem¬ 
perature, and the appropriate time-averaged energy injection rate 
is achieved by adjusting the heating probability. This scheme miti¬ 
gates a common problem in cosmological simulations, namely that 
the energy injected as feedback at each timestep is usually radiated 
away before the gas is able to expand. Moreover, this scheme en¬ 
ables the quantity of energy per feedback event to be specified, ir¬ 
respective of the time-averaged energy injection rate. In the case of 
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feedback associated with star formation, SPH neighbours of newly- 
formed star particles are heated to ATsf = 10^'® K, and the proba¬ 
bility of heating is determined by the fraction of the energy budget 
that is available for feedback, /th- We adopt the convention that 
/th = 1 corresponds to an expectation value of the injected energy 
of 1.736 X 10"^® erg (or 8.73 x ergg“^), the energy 
expected from type II supemovae associated with a Chabrier IMF, 
assuming each SN liberates 10®^ erg. For this energy budget and 
heating temperature, ~ 1.3 fluid elements are expected to be heated 
by each newly-formed star particle (equation|^. 

Four of the simulations presented here, referred to as the ‘cal¬ 
ibrated simulations’ (!|^, adopt different functional forms for /th. 
The efficiency of the subgrid star formation feedback that each sim¬ 
ulation assumes, and the fashion in which this depends upon the 
properties of the local environment, is different in each case. The 
calibrated simulations are realised within a cubic comoving cos¬ 
mological volume of L = 50 cMpc on a side, and each reproduces 
the observed z = 0.1 GSMF to within < 0.3 dex over the range 
10®-® < M* < 10^^ ® Mq (Figure]^. This precision is similar to 
that typically attained by semi-analytic galaxy formation models, 
and is unprecedented for hydrodynamical simulations. The GSMF 
realised in each case is not presented as a prediction, because the 
subgrid feedback parameters were calibrated using the observa¬ 
tions. Nevertheless, the reproduction of the observed GSMF rep¬ 
resents a significant advance. The computational expense of cali¬ 
brating large hydrodynamical simulations precludes an exhaustive 
exploration of the available parameter space, so success was not 
guaranteed, and indeed there was no guarantee that the adopted 
subgrid models could reproduce the observations for any values of 
their parameters. More practically, the availability of hydrodynami¬ 
cal simulations that reproduce the demographics of the local galaxy 
population enables a diverse range of problems to be addressed, 
many of which have previously been limited to investigation with 
semi-analytic models. 

This development adds to recent progress towards the repro¬ 
duction of the observed z = 0.1 GSMF (or, in the case of ‘zoom’ 
simulations, the M* — Mhaio relation) in hydrodynamical simula¬ 
tions (e.g.|Okamoto et al.|2010[ Oppenheimer^e£ay201^ Mun£^ 


al.|2013j Puchwein & Springel|2013[|Vogelsberger et al.|2013[ 


et; 


Hopldns et al.|201^ Vogelsberger et al.|2014a^ . Success in this re- 

gard is typically interpreted as an indication that a particular feed¬ 
back implementation is reasonable. However, the reproduction here 
of the observed z = 0.1 GSMF, with a number of models that yield 
unrealistically compact galaxies, highlights that the simultaneous 
comparison of simulated galaxies with complementary diagnostics, 
such as the observed sizes of galaxies, is necessary to ensure the ac¬ 
curate reproduction of a broad range of galaxy scaling relation|^ 
For example, the FBconst model, which adopts /th = 1 irre¬ 
spective of local conditions, yields a cosmic star formation history 
that is weighted too strongly to early cosmic epochs, resulting in 
the formation of galaxies that are too compact and too quiescent 
at the present epoch. The dominant cause of the shortcomings of 
the FBconst model (and also the FBct and FBZ models) is spu¬ 
rious, numerical radiative losses that artificially suppress the im¬ 
pact of energy feedback. Because (at a fixed resolution and heating 


[Hopkins et al.|jM14) examined the size of galaxies in the FIRE ‘zoom’ 
simulations, concluding they are consistent with observations, whilst the 
Illustris team recently posted a preprint examining their size-mass relation, 
concluding that their ~ L* galaxies are larger than observed by a factor of 
~ 2 jSnyder et al.|20f5|. 


temperature) the cooling time of a hydrodynamic resolution ele¬ 
ment decreases more rapidly with increasing density (oc njj^) than 
does its sound crossing time (oc ), there exists a threshold 

density above which a resolution element radiates away injected 
heat before it can expand ( jPalla Vecchia & Schaye|2012| ). At inter¬ 
mediate resolution, approximately half of all stars in the FBconst 
model form from gas with density above this threshold, resulting 
in the formation of galaxies that are too old, too passive and too 
compact. This demonstrates that it is possible to resolve the over¬ 
cooling problem at the halo scale (i.e. ensuring that haloes form 
galaxies whose masses are consistent with constraints from, for ex¬ 
ample, abundance matching techniques) without addressing the an¬ 
gular momentum problem, and hence still yielding a galaxy popu¬ 
lation that is clearly unrealistic. 

The problem is, partially, a consequence of the inability of 
cosmological simulations to resolve the formation of the first stars. 
The first galaxies whose formation can be captured are necessar¬ 
ily associated with haloes that have not yet been subject to inter¬ 
nal feedback processes, exhibiting unrealistically high gas fractions 
and star formation efficiencies, and potentially triggering a cycle of 
overestimated radiative losses. The FBa and FBZ simulations, for 
which /th increases with decreasing dark matter velocity disper¬ 
sion and metallicity, respectively, successfully interrupt this cycle 
by appealing to /th > 1 for stars forming within nascent galax¬ 
ies. However, this modification alone is insufficient to inhibit the 
onset of numerical overcooling. The comoving stellar mass den¬ 
sity of these simulations rises much more rapidly than observed 
at early cosmic epochs, indicating that star formation is initially 
too efficient, and must decline (too) strongly at later times in order 
to reproduce the observed stellar mass density, and the GSMF, at 
z = 0.1. As in FBconst, galaxies in these simulations exhibit mas¬ 
sive, compact bulges that are a classic signature of overcooling. 

The suppression of numerical losses is crucial not only to en¬ 
sure that the growth of galaxies is appropriately regulated, but also 
to minimise the risk of misinterpretation. If significant numerical 
losses are overlooked, the macroscopic efficiency of feedback (as 
adopted in a particular simulation) will be underestimated. The en¬ 
ergy budget then required for outflows to balance inflows, and so 
to reproduce a realistic galaxy population, will be overestimated. 
Spurious conclusions are then easily drawn. For example, if numer¬ 
ical losses associated with star formation feedback are overlooked, 
one might conclude that the shallow faint end of the GSMF can 
only be reproduced by appealing to non-standard solutions such 
as a top-heavy IMF, gravitational preheating or the adoption of a 
warm dark matter (WDM) cosmogony. In the case of more mas¬ 
sive galaxies, overlooking numerical losses might lead to the in¬ 
ference that very efficient AGN feedback is necessary to regulate 
star formation. These conclusions may of course be reasonable, but 
they should only be inferred from the analysis of simulations with 
demonstrably insignificant numerical radiative losses. 

The ideal solution to the numerical overcooling problem ex¬ 
hibited by the FBconst, FBcr and FBZ simulations is to appeal to 
higher resolution. The critical density for efficient feedback scales 
riH.tc c< rrig ' . Moreover, the improved sampling afforded by 
higher resolution justifies the adoption of a higher heating temper¬ 
ature, which implies a lower heating probability per newly-formed 
star particle. This is an effective means of increasing the critical 
density, since nu.tc cx computational cost of higher 

resolution is currently prohibitive, however, particularly since the 
calibration of subgrid parameters requires that “production” sim¬ 
ulations are preceded by a suite of smaller exploratory runs. Ex- 
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perience from the calibration of the EAGLE simulations indicates 
that volumes of at least L = 25 cMpc are necessary in order to 
adequately sample the population of galaxies of stellar mass corre¬ 
sponding to the knee of the GSME, whose properties are particu¬ 
larly sensitive to the subgrid parameters governing feedback from 
both star formation and AGN. 

Our solution is to adopt a star formation feedback efficiency 
that increases with density at fixed metallicity. This inhibits the 
ISM from reaching the high densities at which our stochastic ther¬ 
mal heating scheme becomes numerically inefficient. This criti¬ 
cal density, ~ 10cm“^, is much greater than the density 

at which a col d interstellar p hase is expected to develop (nn ~ 
0.1 cm“^, e.g. Schaye|2004 1 . We do not include such a phase in 
our simulations, since they lack the resolution required to model 
it accurately. The simulations treat the mutiphase ISM as a single¬ 
phase medium that is subject to a polytropic equation of state; in re¬ 
ality, losses within a cold phase would be significantly higher than 
those expected from this volume-averaged representation, but much 
of the energy injected by feedback is likely to be channelled into 
the hot, diffuse phase of the ISM whose radiative losses will, con¬ 
versely, be much lower. Moreover, the formation of a cold, clumpy 
ISM phase would also lead to more clustered star formation, and 
higher-temperature outflows with relatively low radiative losses. 
Clearly, the simulations cannot model radiative losses accurately 
for gas with nn S> 0.1 cm“^, and the adoption of a density depen¬ 
dence for /th can be justified. 

This scheme enables the EAGLE reference model to be cali¬ 
brated to reproduce the GSME, and the size-mass relation of disc 
galaxies, as observed at z = 0.1. This model therefore resolves 
the overcooling problem at all radii within galaxies, supressing the 
formation of unrealistically massive bulge components and allevi¬ 
ating the angular momentum problem. S15 further demonstrated 
that, with the same calibration, many other observed properties of 
galaxies at z = 0.1 are reproduced, for example specific star for¬ 
mation rates, passive fractions, BH masses, rotation velocities and 
metallicities, in addition to the number density of intergalactic ab¬ 
sorption systems. |Eurlong et al.| ( |20l4l l demonstrated that the EA¬ 
GLE reference model also reproduces the observed evolution of the 
comoving stellar mass density and galaxy specific star formation 
rates. 

Integrating over all feedback events associated with star 
formation, the median subgrid feedback efficiency in the Ref- 
L050N0752 and Ref-L100N1504 simulations is /th = 0.65 and 
0.70, respectively, in spite of this model appealing to fth > 1 when 
stars are bom from high-density gas. The similarity of these values 
to the expected energy budget, i.e. that the median efficiency is of 
order unity, might be considered remarkable. Choosing to calibrate 
the simulations affords the freedom to adopt values much greater or 
much less than unity if such values are required in order to repro¬ 
duce the calibration diagnostics. That such values did not prove to 
be necessary is a non-trivial outcome. It indicates that the radiative 
losses that shape the properties of the galaxy population and their 
environments are established primarily on scales that are “macro¬ 
scopic” in the context of cosmological simulations. 

Having calibrated a model that enables the GSME and sizes of 
galaxies to be reproduced at ^ = 0.1, we have explored the sensi¬ 
tivity of galaxy properties and scaling relations to variation of the 
key parameters of the subgrid routines. In general, the properties of 
the galaxy population are insensitive to the adopted star formation 
threshold. Similarly, when self-regulation is dominated by star for¬ 
mation feedback, the properties of galaxies are largely insensitive 
to the details of the ISM, as was also found by|Schaye et al.|j2010^ 


and |Haas et al.| ( [2013b) . The properties of more massive galaxies, 
whose growth is regulated by AGN feedback, are however sensitive 
to the assumed equation of state of the ISM, because the accretion 
rate onto BHs (and hence the liberation of energy as AGN feed¬ 
back) is inversely proportional to the cube of the local sound speed 
which, for gas subject to a polytropic equation of state, scales as 
Cg OC 


The properties of galaxies are in general sensitive to the value 
of any parameter that affects energetic feedback processes. In corn- 


mon with the findings of many semi-analytic (e.g. 

Benson et al. 

20031|Bower et al.i2006l|Croton et al.|2006l|Somervi 

le et al.|2008 

Bower et al.||2012J and hydrodynamical (e.g. [Haas 

et al.||2013a 

Puchwein & Springel 2013| jVogelsberger et al. 2013 Khandai 

et al.|2014| Vogelsberger et al. 2014al simulations, the properties of 


galaxies with stellar mass < M* at z = 0.1 are largely governed by 
feedback associated with star formation, whilst those of more mas¬ 
sive galaxies are governed by feedback associated with the growth 
of BHs. Feedback regulates the conversion of gas into stars, such 
that, at fixed halo mass, galaxy stellar masses are lower. More ef¬ 
ficient feedback therefore results in the association of galaxies of 
fixed stellar mass with more massive haloes. It also fosters the for¬ 
mation of more extended galaxies, by preferentially ejecting gas 
with the low angular momentum from the ISM. 

Some consequences of more efficient feedback are less intu¬ 
itive. For example, more efficient star formation feedback results 
in higher specific star formation rates at 2 = 0.1. This is because 
galaxies of a fixed stellar mass are hosted by more massive haloes, 
so experience a greater cosmological accretion rate and must in¬ 
crease their SFR to achieve self-regulation. By extension, more ef¬ 
ficient feedback also leads to reduced passive fractions at 2 = 0.1. 
It is also potentially counter-intuitive that more efficient feedback, 
in spite of associating galaxies of fixed stellar mass with more mas¬ 
sive haloes, results in lower central BH masses. This can also be 
explained in terms of self-regulation: more efficient star formation 
feedback enables inflows to be balanced without the need for the 
strong AGN feedback that accompanies rapid BH growth. 

We have not studied simulations in which the subgrid effi¬ 
ciency of AGN feedback, et, has been varied, since, as demon¬ 
strated by |Booth & Schaye] ( |2009||2bl0^ , this only affects the mass 
of BHs. However, effects similar to the variation of /th do accom¬ 
pany the variation of the subgrid viscosity, Cvisc, and the AGN 
heating temperature, ATagn. A higher viscosity enables BHs to 
begin accreting efficiently at a lower mass, boosting the mass of 
BHs at fixed stellar mass, particularly so for the low-mass galax¬ 
ies whose BHs grow primary by viscosity-limited gas accretion 
( |Rosas-Guevara et al.|2013| >. It also results in the formation of more 
extended massive galaxies, since the earlier onset of AGN feedback 
enables the removal of a greater fraction of low angular momentum 
gas from their ISM. A higher heating temperature, despite mak¬ 
ing heating events more intermittent, boosts the efficiency of AGN 
feedback in a fashion similar to increasing fth for star formation 
feedback. Because BH mass is established, to first order, by halo 
mass ( [Booth & Schaye|2010| >, the resulting shift in the halo mass 
associated with a galaxies of fixed stellar mass results in a simi¬ 
lar shift of BH masses, and fosters the formation of more extended 
galaxies. 

Ultimately, the corroboration of our conclusions requires sim¬ 
ulations with sufficient resolution and physical detail to enable ab 
initio prediction of energy and momentum losses within the mul¬ 
tiphase ISM, and the development of such calculations is likely 
to remain beyond our reach in the near future. Whilst progress is 
being made in pursuit of this goal (e.g. [Yirak et al.|2010[ jAluzasj 
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|et al.|[20T2l ), and valuable attempts are underway to bridge the 
vast dynamic range between star-forming complexes and galaxies 


(e.g. 

Joung & Mac Low|20061 Powell et al.|2011 [Hopkins et al.| 

2012 

Agertz et al.|2013 |Creasey et al. 2013[ Hopkins et al.|2014 

Creasey et al. 120151, there are many profitable lines of enquiry to 


pursue with cosmological simulations that appeal to subgrid mod¬ 
els. As was discussed by S15, the chief purpose of such simula¬ 
tions is not to yield ah initio predictions of the masses of galaxies 
and their BHs, but to illuminate the fashion by which astrophysical 
processes shape the evolution of galaxies and their environments. 
Having calibrated simulations to reproduce judiciously chosen ob¬ 
servables, their confrontation with myriad complementary observ¬ 
ables, many more than explored here, represents an efficient and di¬ 
rect means of identifying shortcomings in galaxy formation theory, 
and motivating improved treatments of the relevant astrophysical 
processes. 
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